Elasticity of bilayer lipid membranes from their density correlation function

Abstract

We study the density correlation function (DCF) of DPPC lipid bilayers. We compare Molecular Dynamics (MD) results with theoretical predictions obtained with a mesoscopic description, in terms of the lipid membrane elasticity. One key objective of our work is the quantification of the lipid membrane elasticity directly from the DCF, both for the membrane undulations and local membrane thickness. Our method does not require the definition of instantaneous surfaces or internal variables defining lipid orientations. Building on our previous work, here we focus on the intralayer correlations, i.e. the DCF of lipids residing on the same monolayer, by tracking only the position of the phosphorus atoms in a lipid head group. We demonstrate the relevance of the intralayer two-dimensional (2D) correlations to the total DCF. We further show that all-atom (AA) and coarse grained (CG) lipid forcefields, feature distintively different DCFs. The CG forcefield predicts results in good agreement with the mesoscopic predictions, for the entire wavevector range; the AA forcefield (CHARMM36) predict strong peristaltic fluctuations at long wavevectors $q\gtrsim 0.8$ nm${}^{-1}$, which are absent in the CG lipid model (MARTINI).

GRAPHICAL ABSTRACT

Keywords:

• Bilayer lipid membranes
• capillary waves
• molecular dynamics
• structure factor
• X-ray surface diffraction

1. Introduction

Phospholipid bilayers are quasi-two dimensional (2D) structures that self-assemble spontaneously in water. Lipid membranes play a crucial role in biology, providing a selective permeable layer to cells. Furthermore, lipid membranes have emerged as important artificial soft-matter constructs, which allow performing experiments targeting biophysical relevant questions, by exploiting their flexibility and ability to encapsulate complex fluids. The lipid composition defines the shape and fluctuations of membranes at equilibrium conditions, and under externally applied forces. At length scales much larger than the typical bilayer thickness ($d\approx 5\text{nm}$) the shape of the membrane is usually described as a mathematical surface, $z=\xi (\overrightarrow{x})$ in Monge's representation, with $\overrightarrow{x}\equiv (x,y)$, and projected area $A_o$, on the mean plane of the membrane. The membrane free energy can be defined as a functional of the (assumed) smooth membrane shape, using the Helfrich surface Hamiltonian (HsH) [1], in terms of the mean (κ) and Gaussian (${\kappa }_{\mathrm{G}}$) bending moduli. A spontaneous curvature parameter may be added to describe asymmetric bilayers [2], but here we will restrict ourselves to the simplest case with two lipid monolayers of the same composition. In free standing membranes the surface area self-adjusts to give the tensionless state, with zero surface tension (${\gamma }_o=0$). In experimental situations ${\gamma }_o>0$ may be achieved by applying lateral stress on the membrane perimeter, or by adjusting the osmotic pressure in closed vesicles; in those cases ${\gamma }_o$ must be included in the HsH description as an additional parameter.

The HsH description of a lipid bilayer as a smooth mathematical surface $z=\xi (\overrightarrow{x})$ is accurate for the softer fluctuations modes, i.e. those with low wavevector q in the Fourier representation $\xi (\overrightarrow{x})=\sum _q{\hat{\xi }}_q{\mathrm{e}}^{i\overrightarrow{q}\cdot \overrightarrow{x}}$. These modes feature the largest response to external forces and the largest fluctuations of the bilayer membrane at thermal equilibrium conditions. For the tensionless (${\gamma }_o=0$) state HsH predicts very large mean square amplitude for the Fourier components with low q, as $⟨|{\hat{\xi }}_q{|}^2⟩=\mathrm{kT}/(\kappa q^4{A}_o)$, where T and k are the temperature and Boltzmann constant, respectively.

The analysis of the mean square fluctuations $⟨|{\hat{\xi }}_q{|}^2⟩$ in Molecular Dynamics (MD) simulations allows exploring the relationship between the molecular structure of the lipids and the bending modulus of the bilayer membrane [3–5]. Specifically, the bending modulus can be computed in the tensionless state using $\kappa =\mathrm{kT}/(⟨|{\hat{\xi }}_q{|}^2⟩q^4{A}_o)$ for sufficiently low $|\overrightarrow{q}|$. However, the estimation of κ using this method requires data for very low q values, and hence very large systems, as $q\ge 2\pi /L$, with $L^2=A_o$. Furthermore, proper averaging requires very long MD trajectories, to enable appropriate sampling of the slow fluctuations associated those low-q modes. For higher q, $q^4⟨|{\hat{\xi }}_q{|}^2⟩$ depends on the definition used to construct the mesoscopic surface, $z=\xi (\overrightarrow{x})$, from the analysis of the the instantaneous lipid coordinates obtained throughout the simulated trajectory. Standard definitions might result, when $q\ge 2\pi /L$ is not low enough, in an entanglement of the undulatory (U) mode described by HsH with other internal fluctuation modes of the bilayer membrane, such as protrusions [6].

The traditional approaches to compute κ from MD simulations with relatively low computational effort, rely on the analysis of larger q modes, and that requires to extend the mesoscopic analysis of Helfrich including the fluctuations of the internal structure of the lipid bilayer. In addition to the definition of the membrane shape $\xi (\overrightarrow{x})$, these approaches often require the calculation of other fields, such as the the local bilayer thickness, $d(\overrightarrow{x})$, or the tilt field $\overrightarrow{n}(\overrightarrow{x})$. The latter, describes the local mean orientation of the lipid hydrocarbon chains with respect to the membrane plane [6–9]. The choice of these internal fields and the definitions of the extended Hamiltonian, $\mathcal{H}[\xi (\overrightarrow{x}),d(\overrightarrow{x}),\overrightarrow{n}(\overrightarrow{x}),\dots ]$, is an open question, and a very active area of research. Any mesoscopic description of a molecular system, must achieve a compromise between the internal degrees of freedom needed to keep the accuracy of the Hamiltonian used to fit the MD results, and the numerical method used to extract the elastic constants that enter in $\mathcal{H}$. Usually $\mathcal{H}$ is assumed to be quadratic in the relevant variables, and to converge to the Helfrich Hamiltonian when any field other that $\xi (\overrightarrow{r})$ is maintained to its mean value.

A good extended membrane Hamiltonian allows quantifying κ by analysing small bilayers, i.e. fluctuations with larger q. In addition to bending, these analyses provide valuable information on some internal fluctuations of the lipid membrane. However, the use of the additional mesoscopic fields, introduces both a higher computational cost, and uncertainty on the field definition from the instantaneous molecular configurations. Furthermore, assumptions have to be introduced to combine these fields in the extended Hamiltonian. A comprehensive account of these difficulties, and the disparities in the estimations of the bending modulus of lipid bilayers is given by Ergüder and Deserno [6]. Advancing the discussion below, we note that in simulations of lipid membranes with an all-atom forcefield representation, the role the internal modes is strongly enhanced, with respect to simulations with a coarse grained forcefield for the lipid molecules. Thus, the observed membrane fluctuations in small systems may be very different for different forcefields; whereas in large systems, the undulatory fluctuations are dominant and reflect a bending modulus κ that may be more robust with respect the forcefield choice. This complicates, or even prevents, using a simple extrapolation of the simulation data from small to large sizes, in order to obtain κ; particularly when we try to track the effects of the atomic-scale structure of the lipid molecules.

We recently introduced [10] a deconstruction procedure for the interlayer component of DCF to get κ from the molecular positions, using MD simulations of moderate size lipids bilayers, consisting of a few thousand lipid molecules. Our method circumvents the use of mesoscopic continuous fields $\xi (\overrightarrow{x})$, $d(\overrightarrow{x})$ or $\overrightarrow{n}(\overrightarrow{x})$, and requires only two-site correlations, in terms of the density correlation function (DCF) $G(z_1,z_2,q)$, which is Fourier transformed on the $\overrightarrow{x}$ plane, but maintaining the real-space dependence in the direction normal to the bilayer plane, z. In our analysis of DPPC and POPC phosphatidylcholine lipid bilayers we used only the interlayer component of the DCF between the phosphorus atoms in the lipid head groups. This DCF can be unambiguously computed for any given forcefield representation. The selection of the interlayer DCF, i.e. the correlation between pairs belonging to the opposite lipid layers in the membrane, provides a natural and effective filter for the U fluctuation mode, and it does not require the introduction of additional fields used in the curvature-tilt approaches [6].

In order to advance in the full theoretical description of the simulated $G(z_1,z_2,q)$, and its associated structure factor, we extend in this paper the two-site correlation analysis, to track both the interlayer and the intralayer components of the DCF $G(z_1,z_2,q)$. Again we rely on the tracking of only the phosphorus atoms in the lipid head groups. The idea of the approach presented here, is to extract not only the coupled undulatory mode, with the bending modulus of the lipid bilayer κ, but also to quantitatively characterise the internal fluctuations of the bilayer that may be relevant to analyse, e.g. the impact of proteins embedded in the lipid membrane. The analysis based on the DCF is compared with the mesoscopic description (Section II) in terms of the surfaces $z={\xi }^{\pm }(\overrightarrow{x})$ describing the instantaneous shapes of the ‘upper’ (+) and ‘lower’ (−) lipid monolayers in the bilayer, as calculated with a (well tested) molecular scale definition. To match the DCF and the fluctuating surfaces representation extracted from the same set of MC configurations, we extend the theoretical analysis of Bedeaux and Weeks [11] (BW) to the case of two coupled fluctuating surfaces (Section III). The BW theory defines the contribution of interfacial fluctuations to the DCF $G(z_1,z_2,q)$, using the derivatives of the (one-site) density profile $\rho (z)$ and the mean square values $⟨|{\xi }_q{|}^2⟩$ of the (assumed) Gaussian distribution of the fluctuating interface. In section IV we analyse the contributions to the simulated DCF $G(z_1,z_2,q)$, which are not described in terms of CW fluctuations (double Bedeaux-Weeks formalism, dBW). Specifically, we will focus on the role of the lipid layer 2D compressibility in the intralayer part of $G(z_1,z_2,q)$. Taking into account these non-CW fluctuations (that we refer as ‘correlation background’), in Section V we probe that it is possible to retrieve the simulation values for DCF and the structure factor from its theoretical dBW values. As an extensions of our previous article [10], where only analysed the interlayer component of DCF, we discuss in section VI a deconstruction procedure to obtain the elastic properties of the monolayer from the intralayer DCF by subtracting contributions due to the compressibility of each layer.

Finally, we will discuss the significant differences of $G(z_1,z_2,q)$ for coarse grained (CG) MARTINI [12] and all-atom CHARMM36 [13] forcefields. We will focus on the analysis of internal fluctuations in the bilayer, and demonstrate that both forcefields predict a qualitatively distinct membrane behaviour.

2. Mesoscopic description of the membrane as two coupled fluctuating surfaces

For the analysis presented in this work we employed the MD trajectories used in our previous works of DPPC bilayers [10, 14]. These works included simulations performed with both, all-atom CHARMM36 [13] and CG MARTINI [12] forcefield representations of DPPC molecules in water (see Supplementary Material (SM) for a summary of the simulation details). In the analysis reported below we use only the coordinates ${\overrightarrow{r}}_i^{\alpha }=({\overrightarrow{x}}_i^{\alpha },z_i^{\alpha })$ of the phosphorous atom (i = 1, ··· , $N^{\alpha }$) in each phospholipid head group, for the two paired monolayers ($\alpha =\pm$) with opposite orientations of their hydrophobic tails. The initially equal numbers of phospholipid molecules in each layer, $N^{+}=N^{-}=N/2$, are usually unchanged along the whole MD, since phospholipid flip-flop events are not observed in standard simulations. Nevertheless, if necessary, it is easy to perform some direct test for the head-tail orientation of each phospholipid, to detect any eventual flip-flop event and to relabel that molecule. Any random drift of the membrane along the Z direction is eliminated by taking the origin at the centre of mass, i.e. we take $z_{\mathrm{cm}}\equiv \sum _i(z_i^{+}+z_i^{-})/N=0$ for each configuration, where N is the total number of phospholipids.

The mesoscopic description is used here to compare with our analysis of the membrane elasticity using just the DCF. For each MD configuration we compute the two surfaces $z={\xi }^{+}(\overrightarrow{x})$ and $z={\xi }^{-}(\overrightarrow{x})$, as a smoothed representations of all the ${\overrightarrow{r}}_i^{\alpha }$ in each lipid layer. The undulatory (U) mode of the membrane as a whole is defined with the mean position between the two layers ${\xi }^{\mathrm{U}}(\overrightarrow{x})=({\xi }^{+}(\overrightarrow{x})+{\xi }^{-}(\overrightarrow{x}))/2$ and the peristaltic (P) mode ${\xi }^{\mathrm{P}}(\overrightarrow{x})=({\xi }^{+}(\overrightarrow{x})-{\xi }^{-}(\overrightarrow{x}))/2$ gives the local half-width of the bilayer. In a symmetric bilayer the statistical properties of each monolayer (m) are equal and the functions ${\xi }^{\mathrm{U}}(\overrightarrow{x})$ and ${\xi }^{\mathrm{P}}(\overrightarrow{x})$ become uncorrelated, i.e. $⟨{\xi }^{\mathrm{U}}(\overrightarrow{x}){\xi }^{\mathrm{P}}({\overrightarrow{x}}^{\prime })⟩=0$ for any $\overrightarrow{x}$ and ${\overrightarrow{x}}^{\prime }$, a property that will be very useful for our later analysis.

We refer the reader to our previous works [10, 14] for details of the MD simulations and the mesoscopic analysis to get the monolayers (m) mean square fluctuations $⟨|{\hat{\xi }}_{\overrightarrow{q}}^{\mathrm{m}}{|}^2⟩\equiv (⟨|{\hat{\xi }}_{\overrightarrow{q}}^{+}{|}^2⟩+⟨|{\hat{\xi }}_{\overrightarrow{q}}^{-}{|}^2⟩)/2\approx ⟨|{\hat{\xi }}_{\overrightarrow{q}}^{+}{|}^2⟩\approx ⟨|{\hat{\xi }}_{\overrightarrow{q}}^{-}{|}^2⟩$, and their coupled undulations (CU) $⟨{\hat{\xi }}_{\overrightarrow{q}}^{+}{\hat{\xi }}_{-\overrightarrow{q}}^{-}⟩$. The mean square amplitudes of the normal (U and P) modes are $⟨|{\hat{\xi }}_{\overrightarrow{q}}^{\mathrm{U}}{|}^2⟩=(⟨|{\hat{\xi }}_{\overrightarrow{q}}^{\mathrm{m}}{|}^2⟩+⟨{\hat{\xi }}_{\overrightarrow{q}}^{+}{\hat{\xi }}_{-\overrightarrow{q}}^{-}⟩)/2$ and $⟨|{\hat{\xi }}_{\overrightarrow{q}}^{\mathrm{P}}{|}^2⟩=(⟨|{\hat{\xi }}_{\overrightarrow{q}}^{\mathrm{m}}{|}^2⟩-⟨{\hat{\xi }}_{\overrightarrow{q}}^{+}{\hat{\xi }}_{-\overrightarrow{q}}^{-}⟩)/2$. Figure 1 presents these results obtained with CHARMM36, in terms of the wavevector dependent surface tensions, (1) ${\gamma }^{\mathrm{X}}(q)\equiv \frac{\mathrm{kT}}{⟨|{\hat{\xi }}_q^{\mathrm{X}}{|}^2⟩q^2{A}_o},$(1) for X = m, U and the equivalent with $⟨{\hat{\xi }}_{\overrightarrow{q}}^{+}{\hat{\xi }}_{-\overrightarrow{q}}^{-}⟩$ for CU. The low q regime ${\gamma }^{\mathrm{m}}(q)\approx {\gamma }^{\mathrm{U}}(q)\approx {\gamma }^{\mathrm{CU}}(q)\approx \kappa q^2$ for $q\lesssim 0.35{\text{nm}}^{-1}$, corresponds to the soft modes of the low wavevector undulation for a tensionless membrane, as described in HsH with ${\gamma }_o=0$. Our estimate [10] for the membrane bending modulus is $\kappa /(\mathrm{kT})=26.2\pm 1.5$ [10]. However, the three surface tension functions shown in Figure 1, become qualitatively different for $q\gtrsim 0.5{\text{nm}}^{-1}$, highlighting the relevance of the internal fluctuations of the membrane.

Figure 1. MD simulation results for the variation of the surface tensions with the wavevector q, as defined in (1(1) ${\gamma }^{\mathrm{X}}(q)\equiv \frac{\mathrm{kT}}{⟨|{\hat{\xi }}_q^{\mathrm{X}}{|}^2⟩q^2{A}_o},$(1) ), for the free DPPC membrane modelled with the CHARMM36 forcefield. Green: the coupled undulatory mode ${\gamma }^{\mathrm{CU}}(q)$; red: the peristaltic mode ${\gamma }^{\mathrm{P}}(q)$; blue: the undulatory mode ${\gamma }^{\mathrm{U}}(q)$, and black: the monolayer mode ${\gamma }^{\mathrm{m}}(q)$. The circles correspond to the values of q accessible in the MD box, the lines are guides to the eye.

Within our representation the internal fluctuations emerge from the peristaltic mode, also shown in Figure 1 as ${\gamma }^{\mathrm{P}}(q)\equiv \mathrm{kT}/(⟨|{\hat{\xi }}_q^{\mathrm{P}}{|}^2⟩q^2{A}_o)$. At low q this representation is rather counterintuitive, diverging as $q^{-2}$, since $⟨|{\hat{\xi }}_q^{\mathrm{P}}{|}^2⟩$ has finite q = 0 limit, rather than showing the $q^{-2}$ divergence of any undulatory fluctuation. However, the representation of the P mode as a surface tension function becomes clearly more intuitive for large q, when ${\gamma }^{\mathrm{P}}(q)$ goes smoothly to the same values as the undulatory ${\gamma }^{\mathrm{U}}(q)$, which is about twice that of the monolayers in the high-q regime. This result indicates, since we use only the two surfaces $z={\xi }^{\pm }(\overrightarrow{x})$, that the four functions ${\gamma }^{\mathrm{X}}(q)$ in Figure 1 are linked by the relationships (2) $\begin{array}{rl}& \frac{1}{{\gamma }^{\mathrm{U}}(q)}+\frac{1}{{\gamma }^{\mathrm{P}}(q)}=\frac{1}{{\gamma }^{\mathrm{m}}(q)}\text{and}\\ & \frac{1}{{\gamma }^{\mathrm{U}}(q)}-\frac{1}{{\gamma }^{\mathrm{P}}(q)}=\frac{1}{{\gamma }^{\mathrm{CU}}(q)}.\end{array}$(2) The decoupling of the fluctuations of the two monolayers, $⟨{\hat{\xi }}_{\overrightarrow{q}}^{+}{\hat{\xi }}_{-\overrightarrow{q}}^{-}⟩\approx 0$ at large q, leads to the limits ${\gamma }^{\mathrm{U}}(q)={\gamma }^{\mathrm{P}}(q)=2{\gamma }^{\mathrm{m}}(q)$.

We find qualitative differences between the MD results for DPPC modelled with either all-atom CHARMM36 or CG MARTINI forcefields. In the all-atom forcefield results, see Figure 1, ${\gamma }^{\mathrm{P}}(q)$ crosses over ${\gamma }^{\mathrm{U}}(q)$, with the P mode becoming softer for $q\gtrsim 0.7{\text{nm}}^{-1}$ In contrast, the MARTINI predicts a U mode that is always softer than the P mode (see Figure 5 in the SM). Namely, bending the MARTINI-DPPC bilayer costs less free energy than changing its thickness by the same amplitude and wavelength; but the opposite happens with the CHARMM36 representation of the same lipid molecules, for undulations of less than $\sim 9\text{nm}$ wavelength. In Figure 2 we demonstrate this important difference between all-atom and CG forcefields, by exploring the q-dependence of the CU amplitude $⟨{\hat{\xi }}_{\overrightarrow{q}}^{+}{\hat{\xi }}_{-\overrightarrow{q}}^{-}⟩$. At low q values the results for all-atom and CG representations agree reasonably well predicting similar positive values for the CU amplitude. However, at high-q, the MARTINI forcefield predicts a monotonic decay towards zero, while the all-atom forcefield predicts a negative correlation for $q\gtrsim 0.7{\text{nm}}^{-1}$.

Figure 2. The correlation between the fluctuating amplitudes of the two layers in the membrane, as a function of the wavevector q, for the tensionless (${\gamma }_o=0$) DPPC membrane with two different forcefield representations: The green line shows the result obtained with the MARTINI coarse-grained forcefield, and the black circles the results for the all-atom CHARMM36 forcefield. In order to study the uncertainty of the all-atom results around the minimum, we have analysed the data using subaverages. The dashed line shows the values obtained using only the first 6760 configurations and the dotted line those obtained using only the last 6760 configurations. As we can see the results are well converged The blue squares are results for the tensionless CHARMM36-DPPC membrane obtained with the deconstruction of the interlayer MD $G^{+-}(z_1,z_2,q)$. The red line represents result for the MARTINI-POPC membrane under tension.

In order to characterise the qualitative differences between the MD results for DPPC modelled with either all-atom CHARMM36 or CG MARTINI forcefields, we have analysed the probability distribution of the local width for typical bilayer configurations. Figure 3 shows that for the MARTINI the distribution is a Gaussian with $⟨{\xi }_{\mathrm{P}}⟩\equiv d/2\approx 2.05\pm 0.1\text{nm}$, while for the CHARMM36 we get wider and distinctly asymmetric distribution. The origin of this asymmetry in the membrane width distribution can be easily seen in the snapshot of a representative configuration of CHARMM6 free DPPC, Figure 4. In some regions of the membrane lipid chains remain more or less straight and parallel, implying a large distance between the heads. However, in other areas of the membrane, the flexible CHARMM36 lipid chains show significant deformations with strong folding of the chain, featuring a shorter distance between the heads. Within our description of the membrane, this change in the distance between the lipid heads is expressed as that the peristaltic mode dominating over the undulatory mode at intermediate wave vectors. With MARTINI, the lipid chains are much more rigid, preventing the strong folding of the chains observed with CHARMM36. Although it is logical to argue that the behaviour of a biological membrane will be better described by the more complex forcefield (CHARMM36), the question of which forcefield better captures the behaviour of lipid chains in biological membranes remains open. It is necessary to compare with experimental results of lipid chain behaviour, as well as simulate other membranes and potentials, which exceeds the scope of this work. In this sense the analyses of the deuterium order parameter has been very helpful in previous studies, to assess the accuracy of AA forcefields [15]. However, we will analyse whether these differences between AA and MARTINI descriptions may be directly observed in their DCF, $G(z_1,z_2,q)$.

Figure 3. The dashed red line shows the probability distribution ${\mathcal{P}}_{\mathrm{P}}({\xi }^{\mathrm{P}})$ for the local half-width ${\xi }^{\mathrm{P}}(\overrightarrow{x})$ averaged over all configurations of the free DPPC membranes. ${\xi }^{\mathrm{P}}(\overrightarrow{x})$ is evaluated using the ISM surfaces with $q_u=2.60{\text{nm}}^{-1}$. Top panel: results for the CHARMM36, bottom panel: for the MARTINI. The blue line is the fit with a Gaussian, and (in the top panel) the green line is a fit with two Gaussians.

3. Density profiles and density correlation functions of lipid bilayers

We discuss now the analysis of the same MD trajectories, targeting the density profiles and density correlation function (DCF) of the phosphorous atoms. Here we do not use any mesoscopic representation, such as the surfaces $z={\xi }^{\pm }(\overrightarrow{x})$ discussed above. In our previous work [10] we showed that bending modulus κ can be computed accurately by analysing the interlayer DCF. Here, we explore the use of both, interlayer and intralayer components of the DCF. We investigate their relationship with the internal membrane fluctuations, including the peristaltic mode associated to the membrane thickness, as well as other molecular fluctuation modes, not included in the mesoscopic representation provided by ${\xi }^{\pm }(\overrightarrow{x})$.

Figure 4. Snapshot of a representative configuration of DPPC modeled with the CHARMM36 all-atom forcefield. The snapshot shows two regions in the lipid bilayer corresponding to phosphate-phosphate distances of 48.6 Å(left, yellow, and green lipids) and 38.4 Å(right: orange, cyan, blue, and purple lipids). The phosphorous atoms in the lipid head group are highlighted as red spheres.

Figure 5. TheCHARMM36 DPPC density-density correlation functions from MD simulation with no tensile stress. Left column q = 0.268 nm${}^{-1}$, middle column q = 0.80 nm${}^{-1}$, right column q = 1.34 nm${}^{-1}$. Bottom row: the simulation $G(z_1,z_2;q)$, Equation (4(4) $\begin{array}{rl}G(z_1,z_2,q)& =⟨\sum _{i\ne j=1}^N\delta ({\overrightarrow{r}}_1-{\overrightarrow{r}}_i)\delta ({\overrightarrow{r}}_2-{\overrightarrow{r}}_j){\mathrm{e}}^{i\overrightarrow{q}\cdot {\overrightarrow{x}}_{\mathrm{ij}}}⟩\\ & +\delta (z_1-z_2)\rho (z_1),\end{array}$(4) ) without the delta term. Second row: the double Bedeaux Weeks $G_{\mathrm{dBW}}(z_1,z_2;q)$ to order n = 20 and using an intrinsic profile with $q_u=2.60$ nm${}^{-1}$ (see Equation (8(8) $\begin{array}{rl}& G_{\mathrm{dBW}}({\overrightarrow{x}}_{12},z_1,z_2)\equiv \sum _{\alpha ,\beta =\pm }{G}_{\mathrm{dBW}}^{\mathrm{\alpha \beta }}({\overrightarrow{x}}_{12},z_1,z_2)\\ & \equiv \sum _{\alpha ,\beta =\pm }(⟨{\varrho }_{\xi }^{\alpha }({\overrightarrow{x}}_1,z_1){\varrho }_{\xi }^{\beta }({\overrightarrow{x}}_2,z_2){⟩}_{\xi }-{\rho }^{\alpha }(z_1){\rho }^{\beta }(z_2))\\ & =\sum _{\alpha ,\beta =\pm }\int \int \mathrm{d}{\xi }_1^{\alpha }\mathrm{d}{\xi }_2^{\beta }{\mathcal{P}}^{\mathrm{\alpha \beta }}({\xi }_1^{\alpha },{\xi }_2^{\beta };x_{12})\\ & \times ({\rho }_{\mathrm{Im}}(\alpha (z_1-{\xi }_1^{\alpha })){\rho }_{\mathrm{Im}}(\beta (z_2-{\xi }_2^{\beta }))\\ & -{\rho }^{\alpha }(z_1){\rho }^{\beta }(z_2)\phantom{(\beta (z_2-{\xi }_2^{\beta }))})\end{array}$(8) )). Third row: the sum $G_{\mathrm{dBW}}(z_1,z_2;q)+G_{2\mathrm{Db}}^{++}(z_1,z_2;q)$. Top row: the sum $G_{\mathrm{dBW}}(z_1,z_2;q)+G_{\mathrm{bg}}^{++}(z_1,z_2;q)$.

3.1. Molecular dynamics results

The density profiles for each phospholipid layer can be obtained from the MD averages (3) ${\rho }^{\alpha }(z)=⟨{\varrho }^{\alpha }(\overrightarrow{r})⟩=\frac{1}{A_o}⟨\sum _{i=1}^{N_{\alpha }}\delta (z-z_i^{\alpha })⟩.$(3) where ${\varrho }^{\alpha }(\overrightarrow{r})\equiv \sum _i\delta ({\overrightarrow{r}}_i^{\alpha }-\overrightarrow{r}))$ is the instantaneous density for the positions ${\overrightarrow{r}}_i^{\alpha }$ of the phosphorus atoms. The density profiles (see Figure 3 in SM) feature narrow distributions, ${\rho }^{+}(z)={\rho }^{-}(-z)$, around $z=\pm ⟨d⟩/2$, $⟨d⟩$ being the mean membrane thickness. The total density profile can be obtained from $\rho (z)={\rho }^{+}(z)+{\rho }^{-}(z)$, i.e. by extending the sum in (3(3) ${\rho }^{\alpha }(z)=⟨{\varrho }^{\alpha }(\overrightarrow{r})⟩=\frac{1}{A_o}⟨\sum _{i=1}^{N_{\alpha }}\delta (z-z_i^{\alpha })⟩.$(3) ) over all the phosphorous atoms in the membrane. For not too large membrane sizes the two peaks in $\rho (z)$ are well separated, as shown by the large $\rho (z)\approx 0$ gap centred around z = 0. It is then easy to split ${\rho }^{+}(z)=\rho (z)$ for z>0 and ${\rho }^{-}(z)=\rho (z)$ for z<0, without keeping track of the index $\alpha =\pm$ associated to each phosphorus atom.

The density correlation function (DCF) of the phosphorus atoms at any two points is defined by $G({\overrightarrow{r}}_1,{\overrightarrow{r}}_2)\equiv ⟨\varrho ({\overrightarrow{r}}_1)\varrho ({\overrightarrow{r}}_2)⟩-⟨\varrho ({\overrightarrow{r}}_1)⟩⟨\varrho ({\overrightarrow{r}}_2)⟩$. We compute its 2D Fourier transform in ${\overrightarrow{x}}_{12}={\overrightarrow{x}}_2-{\overrightarrow{x}}_1$ (4) $\begin{array}{rl}G(z_1,z_2,q)& =⟨\sum _{i\ne j=1}^N\delta ({\overrightarrow{r}}_1-{\overrightarrow{r}}_i)\delta ({\overrightarrow{r}}_2-{\overrightarrow{r}}_j){\mathrm{e}}^{i\overrightarrow{q}\cdot {\overrightarrow{x}}_{\mathrm{ij}}}⟩\\ & +\delta (z_1-z_2)\rho (z_1),\end{array}$(4) for all transverse wavevectors $q\equiv |\overrightarrow{q}|\ne 0$ allowed by the MD box size. The product $\rho (z_1)\rho (z_2)$ must be subtracted from (4(4) $\begin{array}{rl}G(z_1,z_2,q)& =⟨\sum _{i\ne j=1}^N\delta ({\overrightarrow{r}}_1-{\overrightarrow{r}}_i)\delta ({\overrightarrow{r}}_2-{\overrightarrow{r}}_j){\mathrm{e}}^{i\overrightarrow{q}\cdot {\overrightarrow{x}}_{\mathrm{ij}}}⟩\\ & +\delta (z_1-z_2)\rho (z_1),\end{array}$(4) ) to get $G(z_1,z_2,0)$, but for our analysis we will only use $q\ne 0$. The last term in (4(4) $\begin{array}{rl}G(z_1,z_2,q)& =⟨\sum _{i\ne j=1}^N\delta ({\overrightarrow{r}}_1-{\overrightarrow{r}}_i)\delta ({\overrightarrow{r}}_2-{\overrightarrow{r}}_j){\mathrm{e}}^{i\overrightarrow{q}\cdot {\overrightarrow{x}}_{\mathrm{ij}}}⟩\\ & +\delta (z_1-z_2)\rho (z_1),\end{array}$(4) ) is the ideal gas self-correlation, independent of q, from the i = j terms that are excluded in the first term sum. The lower panels of Figure 5 shows $G(z_1,z_2,q)$ (without the delta-function self-correlation term) for all-atom CHARMM36 tensionless bilayers, and for three different values of the transverse wavevector q. For systems of moderate size, such as those employed here, $G(z_1,z_2,q)$ features well separated quadrants for positive/negative values of $z_1$ and $z_2$, i.e. the intralayer $G^{\mathrm{\alpha \beta }}(z_1,z_2,q)$ with $\alpha =\beta =\pm$ and interlayer $\alpha \ne \beta$ density correlations, respectively. For larger simulation systems, $L\gtrsim 50\mathrm{nm}$ (see Figure 5 of SM in ref. [10]), we expect some overlap around $z_{1,2}\approx 0$, but still in these cases we may keep track of the ± labels for each lipid molecule to compute the separate averages for each quadrant. The different scale of the colour maps for each q vector (see panel in the same column) reflects the rapid increase of the DCF in a tensionless membrane as the transverse wavevector decreases. The functional dependence of the DCF divergence, $G(z_1,z_2,q)\propto \mathrm{kT}/(\kappa q^4)$, of the thermal undulations at low q, provides the route to obtain the bending modulus of the bilayer membrane [10]. We note that contributions associated the internal modes in the DCF would be more important at large q-vectors, as the contribution of the membrane undulations becomes less important. Figure 6 of the SM shows that the simulation results for $G(z_1,z_2,q)$ using MARTINI closely resembles CHARMM36 data. The only significant difference appears at the highest q value. In the case of CHARMM36, the off-diagonal (interlayer) elements exhibit a change in sign symmetry compared to its lower q values. Conversely, with MARTINI, these elements maintain the same symmetry throughout.

The Fourier transform of $G(z_1,z_2,q)$ with respect to $z_{12}\equiv z_1-z_2$, integrated over $(z_1+z_2)/2$, gives the structure factor: (5) $\begin{array}{rl}S(q,q_z)& =\frac{1}{N}\int \mathrm{d}{z}_1\int \mathrm{d}{z}_{2}G(z_1,z_2,q){\mathrm{e}}^{i{q}_z{z}_{12}}\\ & =1+\frac{1}{N}⟨\sum _{i\ne j=1}^N{\mathrm{e}}^{i(\overrightarrow{q}\cdot {\overrightarrow{x}}_{\mathrm{ij}}+q_z{z}_{\mathrm{ij}})}⟩\end{array}$(5) where we keep the notation $q=|\overrightarrow{q}|$ for the 2D wavevector on the bilayer XY plane. The upper panels in Figure 6 show $S(q,q_z)-1$ (i.e. taking out the ideal-gas self-correlation term in (4(4) $\begin{array}{rl}G(z_1,z_2,q)& =⟨\sum _{i\ne j=1}^N\delta ({\overrightarrow{r}}_1-{\overrightarrow{r}}_i)\delta ({\overrightarrow{r}}_2-{\overrightarrow{r}}_j){\mathrm{e}}^{i\overrightarrow{q}\cdot {\overrightarrow{x}}_{\mathrm{ij}}}⟩\\ & +\delta (z_1-z_2)\rho (z_1),\end{array}$(4) )) for the DCF presented in Figure 5. In our MD simulations the intralayer $S^{++}(q,q_z)=S^{--}(q,q_z)$ and the interlayer $S^{+-}(q,q_z)=S^{-+}(q,q_z)$ components of the structure factor may be separated, with the sum over i, j in (5(5) $\begin{array}{rl}S(q,q_z)& =\frac{1}{N}\int \mathrm{d}{z}_1\int \mathrm{d}{z}_{2}G(z_1,z_2,q){\mathrm{e}}^{i{q}_z{z}_{12}}\\ & =1+\frac{1}{N}⟨\sum _{i\ne j=1}^N{\mathrm{e}}^{i(\overrightarrow{q}\cdot {\overrightarrow{x}}_{\mathrm{ij}}+q_z{z}_{\mathrm{ij}})}⟩\end{array}$(5) ) restricted to the corresponding layers. The interlayer term varies with $q_z$ as a modulated oscillatory function with period $2\pi /⟨d⟩$. The intralayer term is a smooth function of $q_z$, and decays for large $q_z$, when $2\pi /q_z$ becomes small compared with the narrow width of the ${\rho }^{\pm }(z)$ peaks in the density profile. Notice that from $S(q,q_z)$ we cannot recover the full dependence of $G(z_1,z_2,q)$ with $z_1$ and $z_2$. We also note that $G(z_1,z_2,q)$ cannot be accessed experimentally, but the function $S(q,q_z)$ would be closer to a property that might be measured in diffraction experiments. Notwithstanding, in an experimental situation, the simulated $S(q,q_z)$ for the phosphorus atoms would be blurred in the experimental data by a form factor with the scattering section of the molecules (of all the electrons for $X-$rays and all the nuclei for neutrons), as well as the contributions from the water bath. As for the $G(z_1,z_2,q)$ the simulation $S(q,q_z)$ using the MARTINI is very similar to the obtained with the CHARMM36, not presenting notable qualitative differences, see Figures 2 and 7 of the SM.

Figure 6. The results in the two bottom rows of Figure 5 for the DCF of the tensionless CHARMM36-DPPC membrane are presented in terms of the structure factor. Dashed lines: the MD results without the ideal self-correlation, $S(q,q_z)-1$; full lines: the double Bedeaux-Weeks prediction $S_{\mathrm{dBW}}(q,q_z)$, as function of $q_z$. Top panel q = 0.268 nm${}^{-1}$, middle panel q = 0.80 nm${}^{-1}$, and bottom panel q= 1.34 $q=1.34n{m}^{-1}$. Black lines: total structure factor $S(q,q_z)$, blue lines: interlayer component $S^{+-}(q,q_z)=S^{-+}(q,q_z)$ and red lines: intralayer component $S^{++}(q,q_z)=S^{--}(q,q_z)$.

Figure 7. Thenon-CW contribution (background) of the intralayer density correlation function for the same DPPC membrane as in Figure 5. Left column q = 0.80 nm${}^{-1}$, right column q = 1.34 nm${}^{-1}$. $G_{\mathrm{bg}-\mathrm{MD}}^{++}(z_1,z_2,q)$, from the difference between the MD result (4(4) $\begin{array}{rl}G(z_1,z_2,q)& =⟨\sum _{i\ne j=1}^N\delta ({\overrightarrow{r}}_1-{\overrightarrow{r}}_i)\delta ({\overrightarrow{r}}_2-{\overrightarrow{r}}_j){\mathrm{e}}^{i\overrightarrow{q}\cdot {\overrightarrow{x}}_{\mathrm{ij}}}⟩\\ & +\delta (z_1-z_2)\rho (z_1),\end{array}$(4) ) and the dBW prediction (8(8) $\begin{array}{rl}& G_{\mathrm{dBW}}({\overrightarrow{x}}_{12},z_1,z_2)\equiv \sum _{\alpha ,\beta =\pm }{G}_{\mathrm{dBW}}^{\mathrm{\alpha \beta }}({\overrightarrow{x}}_{12},z_1,z_2)\\ & \equiv \sum _{\alpha ,\beta =\pm }(⟨{\varrho }_{\xi }^{\alpha }({\overrightarrow{x}}_1,z_1){\varrho }_{\xi }^{\beta }({\overrightarrow{x}}_2,z_2){⟩}_{\xi }-{\rho }^{\alpha }(z_1){\rho }^{\beta }(z_2))\\ & =\sum _{\alpha ,\beta =\pm }\int \int \mathrm{d}{\xi }_1^{\alpha }\mathrm{d}{\xi }_2^{\beta }{\mathcal{P}}^{\mathrm{\alpha \beta }}({\xi }_1^{\alpha },{\xi }_2^{\beta };x_{12})\\ & \times ({\rho }_{\mathrm{Im}}(\alpha (z_1-{\xi }_1^{\alpha })){\rho }_{\mathrm{Im}}(\beta (z_2-{\xi }_2^{\beta }))\\ & -{\rho }^{\alpha }(z_1){\rho }^{\beta }(z_2)\phantom{(\beta (z_2-{\xi }_2^{\beta }))})\end{array}$(8) ), presented in the middle row, is to be compared with two different approaches to estimate the correlation background (in the second and fourth rows). Top row: The intrinsic correlation result from the ISM analysis of the MD results (14(14) $G_{\mathrm{Im}}(z_1,z_2,q)=⟨\sum _{i,j=1}^N\delta ({\overrightarrow{r}}_1-{\overrightarrow{\tau }}_i)\delta ({\overrightarrow{r}}_2-{\overrightarrow{\tau }}_j){\mathrm{e}}^{i\overrightarrow{q}\cdot {\overrightarrow{x}}_{\mathrm{ij}}}⟩,$(14) ). Bottom row: Intrinsic correlation from a simpler 2D compressibility hypothesis (16(16) $G_{2\mathrm{DI}}(z_1,z_2,q)=h_{2\mathrm{D}}(q)\sum _{\nu =\pm }{\rho }_I^{\nu }(z_1){\rho }_I^{\nu }(z_2).$(16) ). These two intrinsic versions of the missing correlations are convoluted with the monolayer fluctuations to give the predictions for $G_{\mathrm{bg}-\mathrm{MD}}^{++}(z_1,z_2,q)$ in the second (15(15) $\begin{array}{rl}{G}_{\mathrm{bg}}^{++}(z_1,z_2,q)& \approx \int \mathrm{d}{\xi }^{\mathrm{U}}{\mathcal{P}}_{\mathrm{U}}({\xi }^{\mathrm{U}})\int \mathrm{d}{\xi }^{\mathrm{P}}{\mathcal{P}}_{\mathrm{P}}({\xi }^{\mathrm{P}})\\ & \times G_{\mathrm{Im}}(z_1-{\xi }_1^{\pm },z_2-{\xi }_2^{\pm },q),\end{array}$(15) ) and fourth (17(17) $\begin{array}{rl}{G}_{2\mathrm{Db}}^{++}(z_1,z_2,q)& =\int \mathrm{d}{\xi }_m\mathcal{P}({\xi }_m)G_{2\mathrm{DI}}(z_1-\xi ,z_2-\xi ,q)\\ & =h_{2\mathrm{D}}(q)\int \mathrm{d}{\xi }_m\mathcal{P}({\xi }_m)\\ & \times {\rho }_I(z_1-{\xi }_m){\rho }_I(z_2-{\xi }_m)\end{array}$(17) ) rows, respectively.

3.2. Mesoscopic description of the fluctuations effects

A key aspect to extract the elastic properties of the lipid bilayer from ${\rho }_{\alpha }(z)$ and $G^{\mathrm{\alpha \beta }}(z_1,z_2,q)$ is their dependence with the system size. Large membranes develop larger shape fluctuations, which give broader peaks in $\rho (z)$ and broader quadrants in $G(z_1,z_2,q)$ when the density calculations are projected on the plane of the bilayer (normal direction to z in our case). The Capillary Wave Theory (CWT) [16, 17], originally formulated for liquid surfaces, provides the theoretical mesoscopic framework to study this size effects. The theory assumes that the local fluctuations of the position of the interface $z=\xi (\overrightarrow{x})=\sum _{\overrightarrow{q}}{\hat{\xi }}_{\overrightarrow{q}}{\mathrm{e}}^{i\overrightarrow{q}\cdot \overrightarrow{x}}$, have an associated intrinsic density profile ${\rho }_{\mathrm{I}}(z)$ that represents the structure of the liquid surface over a small region, whose size does not depend of the system size. Assuming independent Gaussian fluctuations for the Fourier components ${\hat{\xi }}_q$ and the statistical independence of $\xi (\overrightarrow{x})$ and ${\rho }_{\mathrm{I}}(z)$, a full theoretical framework can be developed in terms of the mean square thermal fluctuations $⟨|{\hat{\xi }}_q{|}^2⟩$ to predict the size and q dependence of $\rho (z)$ and $G(z_1,z_2,q)$. Furthermore, this framework provides a route to extract the elastic properties of the fluctuating surface directly from its DCF.

The reader unfamiliar with the CWT and its practical use in computer simulations through the Intrinsic Sampling Method (ISM) [18] is directed to previous works, both for liquid surfaces [19, 20] and for lipid bilayer membranes [10, 14]. Here we focus on the extension of the theoretical methods to describe the DCF in lipid membranes by considering the mesoscopic fluctuations of the lipid monolayers, described by the two surfaces $z={\xi }^{\pm }(\overrightarrow{x})$. From the ± symmetry of our bilayer membranes we get $⟨{\xi }^{\mathrm{U}}({\overrightarrow{x}}_1){\xi }^{\mathrm{P}}({\overrightarrow{x}}_2)⟩=⟨{\xi }^{\mathrm{U}}⟩⟨{\xi }^{\mathrm{P}}⟩$, for any ${\overrightarrow{x}}_1$ and ${\overrightarrow{x}}_2$. Despite the strong coupling between ${\xi }^{+}({\overrightarrow{x}}_1)$ and ${\xi }^{-}({\overrightarrow{x}}_2)$ for short projected distance $x_{12}=|{\overrightarrow{x}}_1-{\overrightarrow{x}}_2|$, the undulatory (U) and peristaltic (P) fluctuation modes are fully uncorrelated. Therefore, the probability distribution for the local position of the two surfaces ${\xi }^{+}(\overrightarrow{x})$ and ${\xi }^{-}(\overrightarrow{x})$, at the same $\overrightarrow{x}$, may be factorised as $\mathcal{P}({\xi }^{+},{\xi }^{-})={\mathcal{P}}_{\mathrm{U}}({\xi }^{\mathrm{U}}){\mathcal{P}}_{\mathrm{P}}({\xi }^{\mathrm{P}})$, in terms of the probability distribution for each mode.

The ISM [18] applied to bilayer lipid membranes [14], provides an approach to define smooth surfaces $z={\xi }^{\alpha }(\overrightarrow{x})$, which interpolate through the instantaneous positions of the phosphorus atoms in each lipid monolayer. The method uses an upper cutoff parameter $|\overrightarrow{q}|\le q_u$ to set the coarse level of the mesoscopic description. From the functions ${\xi }^{\alpha }(\overrightarrow{x})$, calculated from the analysis of a large set of MD frames, we calculate the probability distributions ${\mathcal{P}}_{\mathrm{U}}({\xi }^{\mathrm{U}})$ and ${\mathcal{P}}_{\mathrm{P}}({\xi }^{\mathrm{P}})$. The CWT assumes that both probabilities are Gaussians, with mean values $⟨{\xi }^{\mathrm{U}}⟩=0$ and $⟨{\xi }^{\mathrm{P}}⟩\equiv d/2$, and mean square fluctuations $⟨|{\xi }^{\mathrm{U}}{|}^2⟩$ and $⟨|{\xi }^{\mathrm{P}}{|}^2⟩$. Within the CWT assumptions, the same set of molecular configurations gives the intrinsic monolayer (Im) density profile ${\rho }_{\mathrm{Im}}(z)=⟨{\varrho }^{\pm }(\overrightarrow{x},\pm z\mp {\xi }^{\pm }(\overrightarrow{x}))⟩$, that may be obtained indistinctly from any of the lipid layers or (to get better statistics) from the average of both for monolayers with the same composition. This ${\rho }_{\mathrm{Im}}(z)$ is a very narrow peak around z = 0, much narrower than ${\rho }^{\pm }(z)$ around $z=\pm d/2$, since ${\rho }_{\mathrm{Im}}(z)$ eliminates from the molecular positions the fluctuations of the mesoscopic surfaces $z={\xi }^{\pm }(\overrightarrow{x})$ (see Figure 3 of the SM). ${\rho }_{\mathrm{Im}}(z)$ reflects our choice for how closely to the positions of the phosphorous atoms we define the intrinsic surfaces ${\xi }^{\alpha }(\overrightarrow{x})$. A perfect interpolation may be achieved with a large $q_u$ (i.e. with enough Fourier components), and then we get a delta-function ${\rho }_{\mathrm{Im}}(z)=n_{2\mathrm{D}}\delta (z)$, with the 2D density of the monolayers. However, the CWT assumptions fail at the molecular scales. To achieve quantitative accuracy of the CWT predictions we must use a lower $q_u$, with some intrinsic width in ${\rho }_{\mathrm{Im}}(z)$ that may not be a fully symmetric function, since ${\rho }_{\mathrm{Im}}(z)$ for z>0 corresponds to fluctuations taking the head groups outwards (towards water) of their mean position, while for z<0 the head groups are inwards (towards the other leaflet).

Following the CWT we use, (6) ${\varrho }_{\xi }(\overrightarrow{x},z)\equiv {\rho }_{\mathrm{Im}}(z-{\xi }^{+}(\overrightarrow{x}))+{\rho }_{\mathrm{Im}}(-z+{\xi }^{-}(\overrightarrow{x})),$(6) as the mesoscopic description of the instantaneous density distribution. Then, we may check if the density profile (3(3) ${\rho }^{\alpha }(z)=⟨{\varrho }^{\alpha }(\overrightarrow{r})⟩=\frac{1}{A_o}⟨\sum _{i=1}^{N_{\alpha }}\delta (z-z_i^{\alpha })⟩.$(3) ), obtained as the statistical average over the atomic positions along the MD simulation, is or not accurately recovered from the average of ${\varrho }_{\xi }(\overrightarrow{x},z)$ over the fluctuations of the surfaces (7) $\begin{array}{rl}{\rho }_{\mathrm{CWT}}(z)& =⟨{\varrho }_{\xi }(\overrightarrow{x},z){⟩}_{\xi }\equiv \int \mathrm{d}{\xi }^{+}d{\xi }^{-}\mathcal{P}({\xi }^{+},{\xi }^{-})\\ & \times ({\rho }_{\mathrm{Im}}(z-{\xi }^{+})+{\rho }_{\mathrm{Im}}(-z+{\xi }^{-}))\\ & =\int \mathrm{d}{\xi }^{\mathrm{U}}{\mathcal{P}}_{\mathrm{U}}({\xi }^{\mathrm{U}})\int \mathrm{d}{\xi }^{\mathrm{P}}{\mathcal{P}}_{\mathrm{P}}({\xi }^{\mathrm{P}})\\ & \times [{\rho }_{\mathrm{Im}}(z-\frac{{\xi }^{\mathrm{U}}+{\xi }^{\mathrm{P}}}{2})\\ & +{\rho }_{\mathrm{Im}}(-z+\frac{{\xi }^{\mathrm{U}}-{\xi }^{\mathrm{P}}}{2})].\end{array}$(7) For $q_u\gtrsim 3.5{\text{nm}}^{-1}$ the mesoscopic description becomes too fine-grained and we observe that this CWT prediction fails (see Figure 3 of the SM). For too low $q_u$ we would need very large systems, since we can analyse only the membrane fluctuation for $|\overrightarrow{q}|\le q_u$. Unless otherwise made explicit, we use here $q_u=2.9{\text{nm}}^{-1}$, that gives a good compromise between the accuracy of the theoretical framework and the computational cost.

To extend the analysis by Bedeaux and Weeks (BW) for the DCF, we define a double BW (dBW) mesoscopic representation of the density correlation in real space (8) $\begin{array}{rl}& G_{\mathrm{dBW}}({\overrightarrow{x}}_{12},z_1,z_2)\equiv \sum _{\alpha ,\beta =\pm }{G}_{\mathrm{dBW}}^{\mathrm{\alpha \beta }}({\overrightarrow{x}}_{12},z_1,z_2)\\ & \equiv \sum _{\alpha ,\beta =\pm }(⟨{\varrho }_{\xi }^{\alpha }({\overrightarrow{x}}_1,z_1){\varrho }_{\xi }^{\beta }({\overrightarrow{x}}_2,z_2){⟩}_{\xi }-{\rho }^{\alpha }(z_1){\rho }^{\beta }(z_2))\\ & =\sum _{\alpha ,\beta =\pm }\int \int \mathrm{d}{\xi }_1^{\alpha }\mathrm{d}{\xi }_2^{\beta }{\mathcal{P}}^{\mathrm{\alpha \beta }}({\xi }_1^{\alpha },{\xi }_2^{\beta };x_{12})\\ & \times ({\rho }_{\mathrm{Im}}(\alpha (z_1-{\xi }_1^{\alpha })){\rho }_{\mathrm{Im}}(\beta (z_2-{\xi }_2^{\beta }))\\ & -{\rho }^{\alpha }(z_1){\rho }^{\beta }(z_2)\phantom{(\beta (z_2-{\xi }_2^{\beta }))})\end{array}$(8) where we use ${\xi }_1^{\pm }\equiv {\xi }^{\pm }({\overrightarrow{x}}_1)$ and ${\xi }_2^{\pm }\equiv {\xi }^{\pm }({\overrightarrow{x}}_2)$ as shorthand to the local height of each monolayer. The joint probability distribution for these pairs of intrinsic surfaces heights, ${\mathcal{P}}^{\mathrm{\alpha \beta }}({\xi }_1^{\alpha },{\xi }_2^{\beta };x_{12})$ depends, on the (in-plane) distance $x_{12}$. From a direct extension of the work of Bedeaux and Weeks [11, 21] for a single fluctuating surface, the Fourier transform of the joint probability distribution may be obtained under the assumption of Gaussian independent distributions for each q, in terms of the height-height correlation function ${\mathcal{S}}^{\mathrm{\alpha \beta }}(x_{12})=⟨{\xi }^{\alpha }(x_1){\xi }^{\beta }(x_2)⟩-⟨{\xi }^{\alpha }(x_1)⟩⟨{\xi }^{\beta }(x_2)⟩$.

The generic form for each quadrant in this mesoscopic dBD description of the DCF is a series (9) $G_{\mathrm{dBW}}^{\mathrm{\alpha \beta }}(z_1,z_2,q)=\sum _{n=1}^{\mathrm{\infty }}(\mathrm{\alpha \beta }{)}^n\frac{{\hat{\mathcal{S}}}_n^{\mathrm{\alpha \beta }}(q)}{n!}\frac{{\mathrm{d}}^n\rho (z_1)}{\mathrm{d}{z}_1^n}\frac{{\mathrm{d}}^n\rho (z_2)}{\mathrm{d}{z}_2^n},$(9) with the dependence on $z_{1,2}$ from the n order derivatives of the density profile, and the dependence with q from the Fourier transform of the the n power of the height-height correlation ${\mathcal{S}}^{\mathrm{\alpha \beta }}(x)$ (10) ${\hat{\mathcal{S}}}_n^{\mathrm{\alpha \beta }}(q)\equiv \int {\mathrm{d}}^2\overrightarrow{x}{[{\mathcal{S}}^{\mathrm{\alpha \beta }}(x)]}^n{\mathrm{e}}^{i\overrightarrow{q}\cdot \overrightarrow{x}}.$(10) The n = 1 term, i.e. the Fourier transform of ${\mathcal{S}}^{\mathrm{\alpha \beta }}(x)$, is given by ${\mathcal{S}}^{\mathrm{\alpha \beta }}(q)\equiv {\mathcal{S}}_1^{\mathrm{\alpha \beta }}(q)=A_o⟨{\xi }_{\overrightarrow{q}}^{\alpha }{\xi }_{-\overrightarrow{q}}^{\beta }⟩$ that may be directly related to the elastic constants energy by the equiparticion of the energy. In particular, for very low q we expect that all the terms $(\alpha ,\beta =\pm )$ would become similar and very large, as (11) ${\mathcal{S}}^{\mathrm{\alpha \beta }}(q)=A_o⟨{\xi }_{\overrightarrow{q}}^{\alpha }{\xi }_{-\overrightarrow{q}}^{\beta }⟩\approx \frac{\mathrm{kT}}{{\gamma }_o{q}^2+\kappa q^4+\cdots }$(11) diverging as $q^{-4}$ for tensionless membranes, and as $q^{-2}$ for membranes under tension.

Obviously, the mesoscopic prediction (8(8) $\begin{array}{rl}& G_{\mathrm{dBW}}({\overrightarrow{x}}_{12},z_1,z_2)\equiv \sum _{\alpha ,\beta =\pm }{G}_{\mathrm{dBW}}^{\mathrm{\alpha \beta }}({\overrightarrow{x}}_{12},z_1,z_2)\\ & \equiv \sum _{\alpha ,\beta =\pm }(⟨{\varrho }_{\xi }^{\alpha }({\overrightarrow{x}}_1,z_1){\varrho }_{\xi }^{\beta }({\overrightarrow{x}}_2,z_2){⟩}_{\xi }-{\rho }^{\alpha }(z_1){\rho }^{\beta }(z_2))\\ & =\sum _{\alpha ,\beta =\pm }\int \int \mathrm{d}{\xi }_1^{\alpha }\mathrm{d}{\xi }_2^{\beta }{\mathcal{P}}^{\mathrm{\alpha \beta }}({\xi }_1^{\alpha },{\xi }_2^{\beta };x_{12})\\ & \times ({\rho }_{\mathrm{Im}}(\alpha (z_1-{\xi }_1^{\alpha })){\rho }_{\mathrm{Im}}(\beta (z_2-{\xi }_2^{\beta }))\\ & -{\rho }^{\alpha }(z_1){\rho }^{\beta }(z_2)\phantom{(\beta (z_2-{\xi }_2^{\beta }))})\end{array}$(8) ) contains only the smooth DCF contributions, proportional to the derivatives of the density profile. The delta-function self-correlation term in (4(4) $\begin{array}{rl}G(z_1,z_2,q)& =⟨\sum _{i\ne j=1}^N\delta ({\overrightarrow{r}}_1-{\overrightarrow{r}}_i)\delta ({\overrightarrow{r}}_2-{\overrightarrow{r}}_j){\mathrm{e}}^{i\overrightarrow{q}\cdot {\overrightarrow{x}}_{\mathrm{ij}}}⟩\\ & +\delta (z_1-z_2)\rho (z_1),\end{array}$(4) ), that adds +1 to the structure factor (5(5) $\begin{array}{rl}S(q,q_z)& =\frac{1}{N}\int \mathrm{d}{z}_1\int \mathrm{d}{z}_{2}G(z_1,z_2,q){\mathrm{e}}^{i{q}_z{z}_{12}}\\ & =1+\frac{1}{N}⟨\sum _{i\ne j=1}^N{\mathrm{e}}^{i(\overrightarrow{q}\cdot {\overrightarrow{x}}_{\mathrm{ij}}+q_z{z}_{\mathrm{ij}})}⟩\end{array}$(5) ), is not included. There are also other sources of correlation between the positions ${\overrightarrow{r}}_i^{\pm }$ and ${\overrightarrow{r}}_j^{\pm }$ of different lipid head groups that are not included in $G_{\mathrm{dBW}}(z_1,z_2,q)$, since they do not correspond to the DCF that arise from the rigid local shift of the intrinsic monolayer density profile assumed to get (8(8) $\begin{array}{rl}& G_{\mathrm{dBW}}({\overrightarrow{x}}_{12},z_1,z_2)\equiv \sum _{\alpha ,\beta =\pm }{G}_{\mathrm{dBW}}^{\mathrm{\alpha \beta }}({\overrightarrow{x}}_{12},z_1,z_2)\\ & \equiv \sum _{\alpha ,\beta =\pm }(⟨{\varrho }_{\xi }^{\alpha }({\overrightarrow{x}}_1,z_1){\varrho }_{\xi }^{\beta }({\overrightarrow{x}}_2,z_2){⟩}_{\xi }-{\rho }^{\alpha }(z_1){\rho }^{\beta }(z_2))\\ & =\sum _{\alpha ,\beta =\pm }\int \int \mathrm{d}{\xi }_1^{\alpha }\mathrm{d}{\xi }_2^{\beta }{\mathcal{P}}^{\mathrm{\alpha \beta }}({\xi }_1^{\alpha },{\xi }_2^{\beta };x_{12})\\ & \times ({\rho }_{\mathrm{Im}}(\alpha (z_1-{\xi }_1^{\alpha })){\rho }_{\mathrm{Im}}(\beta (z_2-{\xi }_2^{\beta }))\\ & -{\rho }^{\alpha }(z_1){\rho }^{\beta }(z_2)\phantom{(\beta (z_2-{\xi }_2^{\beta }))})\end{array}$(8) ).Together with the trivial self-correlation term, they give a correlation background $G-G_{\mathrm{dBW}}$ to be discussed in the next subsection.

In Figure 5 the (bottom row) MD results for $G(z_1,z_2,q)$ may be compared with the dBW prediction given in the second from bottom row, when we plug in (9(9) $G_{\mathrm{dBW}}^{\mathrm{\alpha \beta }}(z_1,z_2,q)=\sum _{n=1}^{\mathrm{\infty }}(\mathrm{\alpha \beta }{)}^n\frac{{\hat{\mathcal{S}}}_n^{\mathrm{\alpha \beta }}(q)}{n!}\frac{{\mathrm{d}}^n\rho (z_1)}{\mathrm{d}{z}_1^n}\frac{{\mathrm{d}}^n\rho (z_2)}{\mathrm{d}{z}_2^n},$(9) ) the ISM results for ${\hat{\mathcal{S}}}_n^{\mathrm{\alpha \beta }}(q)$. For the lowest wavevector ($q=0.268{\text{nm}}^{-1}$, on the left column) the visual agreement is quite good. The strong correlations observed in the MD results are clearly the consequence of the large fluctuations of the membrane shape, that are well described by the dBW assumptions. In the middle column (q = 0.80 nm${}^{-1}$) the visual agreement between MD and dBW is still good, including the clear difference between the intralayer (same sign of $z_1$ and $z_2$) and the interlayer ($z_1$ and $z_2$ with opposite sign) quadrants, which reflects the strong contribution to the DCF from the local fluctuations of the membrane width. Still, with a closer look, we may observe that the dBW predictions overestimates the positive values of the MD results $G(z_1,z_2,q)$ for $z_1\approx z_2$ (intralayer). In the right column, for q = 1.34 nm${}^{-1}$ , the difference between the intra and the (much fainter) interlayer quadrants increases. Although the dBW prediction still captures the dominant behaviour, the close comparison with the MD results shows that dBW underestimates the weak interlayer DCF, while it clearly overestimates the stronger intralayer one. On the contrary, for the MARTINI results, at these high q values, the mesoscopic values of $G_{\mathrm{dBW}}^{\pm ,\mp }$ continue to be very similar to the simulation values (see Figure 6 in the SM).

To better quantify the accuracy of dBW, we compare the MD structure factor (5(5) $\begin{array}{rl}S(q,q_z)& =\frac{1}{N}\int \mathrm{d}{z}_1\int \mathrm{d}{z}_{2}G(z_1,z_2,q){\mathrm{e}}^{i{q}_z{z}_{12}}\\ & =1+\frac{1}{N}⟨\sum _{i\ne j=1}^N{\mathrm{e}}^{i(\overrightarrow{q}\cdot {\overrightarrow{x}}_{\mathrm{ij}}+q_z{z}_{\mathrm{ij}})}⟩\end{array}$(5) ), as $S(q,q_z)-1$, to remove the ideal gas self-correlation, with its mesoscopic dBW form (8(8) $\begin{array}{rl}& G_{\mathrm{dBW}}({\overrightarrow{x}}_{12},z_1,z_2)\equiv \sum _{\alpha ,\beta =\pm }{G}_{\mathrm{dBW}}^{\mathrm{\alpha \beta }}({\overrightarrow{x}}_{12},z_1,z_2)\\ & \equiv \sum _{\alpha ,\beta =\pm }(⟨{\varrho }_{\xi }^{\alpha }({\overrightarrow{x}}_1,z_1){\varrho }_{\xi }^{\beta }({\overrightarrow{x}}_2,z_2){⟩}_{\xi }-{\rho }^{\alpha }(z_1){\rho }^{\beta }(z_2))\\ & =\sum _{\alpha ,\beta =\pm }\int \int \mathrm{d}{\xi }_1^{\alpha }\mathrm{d}{\xi }_2^{\beta }{\mathcal{P}}^{\mathrm{\alpha \beta }}({\xi }_1^{\alpha },{\xi }_2^{\beta };x_{12})\\ & \times ({\rho }_{\mathrm{Im}}(\alpha (z_1-{\xi }_1^{\alpha })){\rho }_{\mathrm{Im}}(\beta (z_2-{\xi }_2^{\beta }))\\ & -{\rho }^{\alpha }(z_1){\rho }^{\beta }(z_2)\phantom{(\beta (z_2-{\xi }_2^{\beta }))})\end{array}$(8) ) we get $\begin{array}{rl}{S}_{\mathrm{dBW}}^{\mathrm{\alpha \beta }}(q,q_z)& ={\mathrm{\Phi }}^{\alpha }(q_z){\mathrm{\Phi }}^{\beta }(-q_z)\sum _{n=1}^{\mathrm{\infty }}(\mathrm{\alpha \beta }{q}_z^2{)}^n\frac{{\hat{\mathcal{S}}}_n^{\mathrm{\alpha \beta }}(q)}{n!}\\ & ={\mathrm{\Phi }}^{\alpha }(q_z){\mathrm{\Phi }}^{\beta }(-q_z)\int {\mathrm{d}}^2\overrightarrow{x}\\ & \times ({\mathrm{e}}^{\mathrm{\alpha \beta }{q}_z^2{\mathcal{S}}^{\mathrm{\alpha \beta }}(x)}-1){\mathrm{e}}^{i\overrightarrow{q}\cdot \overrightarrow{x}}.\end{array}$The MD results for the surface structure factor are well approximated by the CWT representation of the mean monolayer profile as a Gaussian convolution with square width ${\mathrm{\Delta }}_m$ of an (very narrow) intrinsic monolayer profile with 2D density $n_{2\mathrm{D}}$; that gives: (12) ${\mathrm{\Phi }}^{\pm }(q_z)=\int \mathrm{d}z\frac{\mathrm{d}{\rho }^{\pm }(z)}{\mathrm{d}z}{\mathrm{e}}^{i{q}_{z}z}\approx -i{n}_{2\mathrm{D}}{\mathrm{e}}^{-\frac{{\mathrm{\Delta }}_m}{2}{q}_z^2}{\mathrm{e}}^{\pm i\frac{⟨d⟩}{2}{q}_z}.$(12) Therefore $|{\mathrm{\Phi }}^{\pm }(q_z){|}^2\approx q_z^2{n}_{2\mathrm{D}}^2{\mathrm{e}}^{-{\mathrm{\Delta }}_{\mathrm{m}}{q}_z^2}$ for the intralayer terms, with a smooth shape for $S_{\mathrm{dBW}}^{\mathrm{\alpha \alpha }}(q,q_z)$ vanishing both at low and high $q_z$. For the interlayer components the sum of the two (off-diagonal) terms gives (13) $\begin{array}{rl}& {\mathrm{\Phi }}^{+}(q_z){\mathrm{\Phi }}^{-}(-q_z)+{\mathrm{\Phi }}^{-}(q_z){\mathrm{\Phi }}^{+}(-q_z)\\ & \approx (|{\mathrm{\Phi }}^{+}(q_z){|}^2+|{\mathrm{\Phi }}^{-}(q_z){|}^2)\cos (d{q}_z),\end{array}$(13) with an oscillatory shape convoluted by the same smooth shape of the intralayer components.

The top panel in Figure 6 shows that for our lowest wavevector (q = 0.268 nm${}^{-1}$) $S_{\mathrm{dBW}}(q,q_z)$ is large and nearly indistinguishable from the MD results. The lower scale for q = 0.80 nm${}^{-1}$ makes visible some discrepancy for the intralayer contribution at low $q_z$, and with the still lower scale for q = 1.34 nm${}^{-1}$ the discrepancy becomes obvious and also extended to higher $q_z$. The qualitative failure of (12(12) ${\mathrm{\Phi }}^{\pm }(q_z)=\int \mathrm{d}z\frac{\mathrm{d}{\rho }^{\pm }(z)}{\mathrm{d}z}{\mathrm{e}}^{i{q}_{z}z}\approx -i{n}_{2\mathrm{D}}{\mathrm{e}}^{-\frac{{\mathrm{\Delta }}_m}{2}{q}_z^2}{\mathrm{e}}^{\pm i\frac{⟨d⟩}{2}{q}_z}.$(12) ) comes in the prediction $S_{\mathrm{dBW}}^{\mathrm{\alpha \alpha }}(q,q_z)=0$ at $q_z=0$, while the MD results for $S^{\mathrm{\alpha \alpha }}(q,q_z)-1\approx -1$, at low $q_z$. Notice that $q_z=0$ corresponds to the integral across the z direction, so that $S^{\mathrm{\alpha \alpha }}(q,0)$ represents the 2D structure factor of the lipid heads in the same monolayer if projected on the XY plane. That contribution to the DCF is obviously out of the dBW description that considers only the effects of the rigid shift of the intrinsic monolayer density profiles ${\rho }_{\mathrm{Im}}(z)$, without any in-plane fluctuation of the intrinsic density. The MD interlayer contribution is $S^{+-}(q,q_z)\approx 0$ for low $q_z$ (the self-correlation does not appear, so that we do not have to take out the $-1$ term) and it compares much better with the dBW result at q = 0.80 nm${}^{-1}$ , although with some phase shift between the MD and the dBW results.

Our previous proposal to extract the bending modulus of a lipid bilayer membrane is based on the identification of the MD result for the interlayer DFC with the dBW prediction (9(9) $G_{\mathrm{dBW}}^{\mathrm{\alpha \beta }}(z_1,z_2,q)=\sum _{n=1}^{\mathrm{\infty }}(\mathrm{\alpha \beta }{)}^n\frac{{\hat{\mathcal{S}}}_n^{\mathrm{\alpha \beta }}(q)}{n!}\frac{{\mathrm{d}}^n\rho (z_1)}{\mathrm{d}{z}_1^n}\frac{{\mathrm{d}}^n\rho (z_2)}{\mathrm{d}{z}_2^n},$(9) ). Assuming that, for low and middle values of q, we may take $G_{\mathrm{dBW}}^{+-}(z_1,z_2,q)\approx G^{+-}(z_1,z_2,q)$, and together with a Gaussian approximation for the MD density profiles ${\rho }^{+}(z)$ and ${\rho }^{-}(z)$, the dBW series may be deconstructed to get the functions ${\hat{\mathcal{S}}}_n^{+-}(q)$, and the elastic constants of the membrane are obtained from the n = 1 term as in (11(11) ${\mathcal{S}}^{\mathrm{\alpha \beta }}(q)=A_o⟨{\xi }_{\overrightarrow{q}}^{\alpha }{\xi }_{-\overrightarrow{q}}^{\beta }⟩\approx \frac{\mathrm{kT}}{{\gamma }_o{q}^2+\kappa q^4+\cdots }$(11) ). Nevertheless, for the all-atom forcefield of DPPC, we find that the discrepancy between MD and dBW becomes much stronger for larger wavevectors. At q = 1.34 nm${}^{-1}$ (on the bottom panel of Figure 6) $S_{\mathrm{dBW}}^{+-}(q,q_z)$ is much smaller than the MD $S^{+-}(q,q_z)$. That discrepancy corresponds to the fainter colour map of the dBW results in the interlayer quadrants of the right hand column in Figure 5. Notice that, over that region of high q and for the CHARMM36 representation of the DPPC lipids, we are getting the negative correlation between the two lipid layers shown in Figure 2. That is in qualitative contrast with the results obtained with the CG MARTINI forcefield at high q where the correlation between layers keep positive and the $S_{\mathrm{dBW}}^{+-}(q,q_z)$ is very similar to the MD $S^{+-}(q,q_z)$ (see Figure 7 in SM)

4. Analysis of the correlation background

The correlation background $G_{\mathrm{bg}-\mathrm{MD}}(z_1,z_2,q)=G(z_1,z_2,q)-G_{\mathrm{dBW}}(z_1,z_2,q)$ gives the DCF that, as shown in the MD simulations, gives the difference between the DCF obtained directly in the MD simulations and the prediction of the mesoscopic theory. That background becomes relatively more important for larger q and it is rather different for intralayer than for interlayer DCF.

4.1. The intralayer correlation background

We consider first the DCF within a lipid monolayer and its background $G_{\mathrm{bg}-\mathrm{MD}}^{++}(z_1,z_2,q)$ (see the middle row of Figure 7) for the two higher values of q shown in Figure 5. Notice that, in contrast with the need of very different colour-scales for each q used in Figure 5, $G_{\mathrm{bg}-\mathrm{MD}}^{++}(z_1,z_2,q)$ depends weakly on q. The elongated shape along the $z_1=z_2$ diagonal, reflects the width of the monolayer density profile ${\rho }^{+}(z)={\rho }^{-}(-z)$, induced by the fluctuations of the lipid layer, and the shape would become even more elongated for membranes with larger size. The narrower span along the other diagonal illustrates the relevant correlations responsible for the difference between the MD result and the dBW mesoscopic prediction.

4.1.1. MD-ISM intrinsic intralayer correlation background

Within the CWT formalism, the background correlations are the fluctuations existing in the intrinsic structure of the membrane, i.e. how the molecular fluctuations would be observed if we virtually translate all the molecules to give flat $z={\xi }^{\pm }(\overrightarrow{x})$ ISM intrinsic surfaces. These can be evaluated following the same lines used to get the intrinsic density profiles of each monolayer in our ISM-MD analysis. We define for each phosphorus atom the instantaneous intrinsic position, ${\overrightarrow{\tau }}_i=({\overrightarrow{x}}_i,z_i-{\xi }^{\pm }({\overrightarrow{x}}_i))$ shifting the position along with the intrinsic surface ${\xi }^{\pm }(\overrightarrow{x})$ for the layer at which the molecule belongs. The MD average for the two-particle intrinsic positions gives (14) $G_{\mathrm{Im}}(z_1,z_2,q)=⟨\sum _{i,j=1}^N\delta ({\overrightarrow{r}}_1-{\overrightarrow{\tau }}_i)\delta ({\overrightarrow{r}}_2-{\overrightarrow{\tau }}_j){\mathrm{e}}^{i\overrightarrow{q}\cdot {\overrightarrow{x}}_{\mathrm{ij}}}⟩,$(14) for any $q\ne 0$. The results in the top row of Figure 7, show very narrow functions on the $(z_1,z_2)$ plane, since we have eliminated the (size dependent) broadening produced by the fluctuations of the monolayer surfaces. Under the hypothesis of local transverse dependence, which was found to be accurate for liquid interfaces [22], we may estimate the intralayer correlation background as (15) $\begin{array}{rl}{G}_{\mathrm{bg}}^{++}(z_1,z_2,q)& \approx \int \mathrm{d}{\xi }^{\mathrm{U}}{\mathcal{P}}_{\mathrm{U}}({\xi }^{\mathrm{U}})\int \mathrm{d}{\xi }^{\mathrm{P}}{\mathcal{P}}_{\mathrm{P}}({\xi }^{\mathrm{P}})\\ & \times G_{\mathrm{Im}}(z_1-{\xi }_1^{\pm },z_2-{\xi }_2^{\pm },q),\end{array}$(15) with the signs in ${\xi }_{1,2}^{\pm }=({\xi }^{\mathrm{U}}\pm {\xi }^{\mathrm{p}})/2$ taken as those of $z_{1,2}$. The results, presented in the second row of Figure 7 show the broadening of the narrow intrinsic structure through its convolution with the fluctuations of the monolayer. Compared with the difference between the MD and the dBW results (third row of the same figure), this estimation seems to give a good description of $G-G_{\mathrm{dBW}}$ along the main diagonal axis ($z_1=z_2$), but it clearly fails to reproduce the form of correlation background along the other diagonal axis ($z_1-d/2=d/2-z_2$). The problem seems to be associated to the strong anysotropic structure of $G_{\mathrm{Im}}(z_1,z_2,q)$, with negative values along the $z_1\approx z_2$ diagonal, but with positive values along the other diagonal which cannot be dispersed by the convolution (15(15) $\begin{array}{rl}{G}_{\mathrm{bg}}^{++}(z_1,z_2,q)& \approx \int \mathrm{d}{\xi }^{\mathrm{U}}{\mathcal{P}}_{\mathrm{U}}({\xi }^{\mathrm{U}})\int \mathrm{d}{\xi }^{\mathrm{P}}{\mathcal{P}}_{\mathrm{P}}({\xi }^{\mathrm{P}})\\ & \times G_{\mathrm{Im}}(z_1-{\xi }_1^{\pm },z_2-{\xi }_2^{\pm },q),\end{array}$(15) ). The most likely interpretation is that we are reaching the limits for the simplest assumption $\xi ({\overrightarrow{x}}_1)\approx \xi ({\overrightarrow{x}}_2)$, behind the convolution (15(15) $\begin{array}{rl}{G}_{\mathrm{bg}}^{++}(z_1,z_2,q)& \approx \int \mathrm{d}{\xi }^{\mathrm{U}}{\mathcal{P}}_{\mathrm{U}}({\xi }^{\mathrm{U}})\int \mathrm{d}{\xi }^{\mathrm{P}}{\mathcal{P}}_{\mathrm{P}}({\xi }^{\mathrm{P}})\\ & \times G_{\mathrm{Im}}(z_1-{\xi }_1^{\pm },z_2-{\xi }_2^{\pm },q),\end{array}$(15) ), and a consistent evaluation of the background from MD-ISM simulations becomes impracticable beyond that approximation.

4.1.2. Intralayer correlation background as the compressibility of lipid monolayers

The two lower rows in Figure 7 correspond to a simpler assumption for the intralayer density correlations not included in the dBW mesoscopic description. The elliptic negative form of the results for $G-G_{\mathrm{dBW}}$ in the $(z_1,z_2)$ map (third row in Figure 7) could be achieved if the convolution (15(15) $\begin{array}{rl}{G}_{\mathrm{bg}}^{++}(z_1,z_2,q)& \approx \int \mathrm{d}{\xi }^{\mathrm{U}}{\mathcal{P}}_{\mathrm{U}}({\xi }^{\mathrm{U}})\int \mathrm{d}{\xi }^{\mathrm{P}}{\mathcal{P}}_{\mathrm{P}}({\xi }^{\mathrm{P}})\\ & \times G_{\mathrm{Im}}(z_1-{\xi }_1^{\pm },z_2-{\xi }_2^{\pm },q),\end{array}$(15) ) were applied to a simple negative shape of the intrinsic background for the monolayer, rather than to the structured form $G_{\mathrm{Im}}$ given by (14(14) $G_{\mathrm{Im}}(z_1,z_2,q)=⟨\sum _{i,j=1}^N\delta ({\overrightarrow{r}}_1-{\overrightarrow{\tau }}_i)\delta ({\overrightarrow{r}}_2-{\overrightarrow{\tau }}_j){\mathrm{e}}^{i\overrightarrow{q}\cdot {\overrightarrow{x}}_{\mathrm{ij}}}⟩,$(14) ). That is the effect expected from the in-surface compressibility fluctuations of the lipid layers; i.e. the fluctuations in the 2D density of lipid heads within the surfaces $z={\xi }^{\pm }(\overrightarrow{x})$, which are not accounted by fluctuations in the shape of these surfaces. The very low values of the interlayer contribution to $G-G_{\mathrm{dBW}}$ would reflect the weak coupling between the two layers through their flexible hydrocarbon tails, while for the intralayer correlation the steric molecular repulsions should give a strong signal.

From MD simulations we have obtained the (projected) 2D structure factor $S^{2\mathrm{D}}(q)$ and the Fourier transformed total correlation function $h_{2\mathrm{D}}(q)=(S_{2\mathrm{D}}(q)-1)/n_{2\mathrm{D}}$, with the mean density $n_{2\mathrm{D}}$ of the monolayer. Then we approximate the intralayer intrinsic correlation background as (16) $G_{2\mathrm{DI}}(z_1,z_2,q)=h_{2\mathrm{D}}(q)\sum _{\nu =\pm }{\rho }_I^{\nu }(z_1){\rho }_I^{\nu }(z_2).$(16) As in Equation (15(15) $\begin{array}{rl}{G}_{\mathrm{bg}}^{++}(z_1,z_2,q)& \approx \int \mathrm{d}{\xi }^{\mathrm{U}}{\mathcal{P}}_{\mathrm{U}}({\xi }^{\mathrm{U}})\int \mathrm{d}{\xi }^{\mathrm{P}}{\mathcal{P}}_{\mathrm{P}}({\xi }^{\mathrm{P}})\\ & \times G_{\mathrm{Im}}(z_1-{\xi }_1^{\pm },z_2-{\xi }_2^{\pm },q),\end{array}$(15) ), the correlation background broadened by the fluctuations of the monolayer is (17) $\begin{array}{rl}{G}_{2\mathrm{Db}}^{++}(z_1,z_2,q)& =\int \mathrm{d}{\xi }_m\mathcal{P}({\xi }_m)G_{2\mathrm{DI}}(z_1-\xi ,z_2-\xi ,q)\\ & =h_{2\mathrm{D}}(q)\int \mathrm{d}{\xi }_m\mathcal{P}({\xi }_m)\\ & \times {\rho }_I(z_1-{\xi }_m){\rho }_I(z_2-{\xi }_m)\end{array}$(17) where $\mathcal{P}({\xi }_m)$ is the probability of the monolayer fluctuations. The results in the two bottom rows in Figure 7 show that the required shape of $G-G_{\mathrm{dBW}}$ may indeed be obtained from a 2D- structure factor background, under a reasonable and simple hypothesis. As we can see in Figure 8 of the SM, this agreement is even much better when simulating membranes with the MARTINI representation.

Figure 8. The structure factor $S(q,q_z)$ of the free All-Atoms DPPC membrane for at q = 0.80 nm${}^{-1}$ (top panel), and q = 1.34 nm${}^{-1}$ (bottom panel). Black lines: direct MD results, dashed lines total structure factor and full line the intralayer contribution. Red full line: Bedeaux-Weeks (dBW) theoretical prediction to $n_{\mathrm{BW}}=20$ for the intralayer contribution, obtained by Fourier transform of $G_{\mathrm{dBW}}^{++}(z_1,z_2,q)$ (see Equation (12(12) ${\mathrm{\Phi }}^{\pm }(q_z)=\int \mathrm{d}z\frac{\mathrm{d}{\rho }^{\pm }(z)}{\mathrm{d}z}{\mathrm{e}}^{i{q}_{z}z}\approx -i{n}_{2\mathrm{D}}{\mathrm{e}}^{-\frac{{\mathrm{\Delta }}_m}{2}{q}_z^2}{\mathrm{e}}^{\pm i\frac{⟨d⟩}{2}{q}_z}.$(12) )). The orange and green lines are $S_{\mathrm{dBW}}^{++}(q,q_z)+S_{\mathrm{bg}}^{++}(q,q_z)$. Orange lines: $S_{\mathrm{bg}}^{++}(q,q_z)$ from Fourier transform of $G_{\mathrm{bg}}^{++}(z_1,z_2,q)$ and green lines of $G_{2\mathrm{Db}}^{++}(z_1,z_2,q)$. The full line: $q_u=1.88$ nm${}^{-1}$ and dashed line: $q_u=2.95$ nm${}^{-1}$.

5. The interlayer correlation background

For the interlayer component the correlation background $G_{\mathrm{bg}-\mathrm{MD}}^{+-}(z_1,z_2,q)=G^{+-}(z_1,z_2,q)-G_{\mathrm{dBW}}^{+-}(z_1,z_2,q)$ is very small up to $q\lesssim 0.6{\text{nm}}^{-1}$. This fact was used in our previous work [10] to calculate the bending modulus from the assumption that the MD result for $G^{+-}(z_1,z_2,q)$ may be directly interpreted as $G_{\mathrm{dBW}}^{+-}(z_1,z_2,q)$. However, for higher q = 1.34 nm${}^{-1}$ , we observe that, being the interlayer component weaker than the intralayer one in the MD results, it is still much weaker in the mesoscopic dBW prediction, and it cannot be improved by our estimations for $G_{\mathrm{bg}}$ (see in Figure 5 the right column of the third and top row). That observation is apecific to all-atoms forcefield (as shown in Figure 4), since these forcefields allow occasional strong folding of the lipid chains. If we use the coarser MARTINI for the same lipid molecule we get that both the MD results and the mesoscopic dBW estimation give similar interlayer DCF at $q\gtrsim$ 1 nm${}^{-1}$ (see Figure 6 in SM). This difference between the two (all-atom and coarser) representations of the same lipid molecule is clearly associated to their different results in Figure 2. The smooth mononotonic decay of $⟨{\hat{\xi }}_q^{+}{\hat{\xi }}_{-q}^{-}⟩$ for increasing q observed with MARTINI is in qualitative contrast with the change of sign for this correlation observed with the all- atom CHARMM36.

The fact that the observation of negative correlations depends on the forcefield used, highlights the need for experimental verification of whether such correlation do indeed exist. Furthermore our results motivate additional computations with other forcefields. As we discuss below, if this aspect of the DCF, observed with the all-atoms representation and missed by MARTINI, is a true characteristic of bilayer lipid membranes, it would have physical significance for the elasticity of these systems at nanometric scale.

5.1. DFC and structure factor predicted from Bedeaux–Weeks plus correlation background

The two upper rows in Figure 5 present the predictions for the full DCF (intralayer and interlayer quadrants) adding to the mesoscopic $G_{\mathrm{dBW}}$ prediction either one or the other of the above theoretical methods to estimate the correlation background, to be compared with the (bottom row) MD results. For $q\le 0.80$ nm${}^{-1}$ the two methods give similar contributions to the intralayer background, while their contribution to the interlayer correlation is negligible. At the right hand column in the same figure (q = 1.34 nm${}^{-1}$ ), we observe that the intra/interlayer asymmetry is much stronger, which makes much clearer the need of a correlation background beyond the double BW series. The estimation of the background based on the 2D structure factor (included in the second row of Figure 7, $G_{2\mathrm{Db}}^{++}$) is clearly more accurate than the estimation from the monolayer intrinsic correlation (included in the four row of Figure 7, $G_{\mathrm{bg}}^{++}$) under the assumption $⟨{\xi }^{+}({\overrightarrow{x}}_1)⟩=⟨{\xi }^{+}({\overrightarrow{x}}_2)⟩$ used in (15(15) $\begin{array}{rl}{G}_{\mathrm{bg}}^{++}(z_1,z_2,q)& \approx \int \mathrm{d}{\xi }^{\mathrm{U}}{\mathcal{P}}_{\mathrm{U}}({\xi }^{\mathrm{U}})\int \mathrm{d}{\xi }^{\mathrm{P}}{\mathcal{P}}_{\mathrm{P}}({\xi }^{\mathrm{P}})\\ & \times G_{\mathrm{Im}}(z_1-{\xi }_1^{\pm },z_2-{\xi }_2^{\pm },q),\end{array}$(15) ). These results confirm that, beside the U and P modes for the fluctuations in the shape of the lipid layers, the next relevant contribution to the density correlations in a bilayer lipid membrane is the 2D compressibility fluctuations.

In Figure 6 we compared the MD results with the mesoscopic dBW prediction for the structure factor $S(q,q_z)$. In Figure 8 we add the estimations from the correlation background that provide a test of their accuracy, better than the colour maps for $G(z_1,z_2,q)$ in Figure 7. The results are presented only for q = 0.80 and 1.34 nm${}^{-1}$ , since for the lowest wavevector q = 0.268 nm${}^{-1}$ the mesoscopic dBW prediction presented in Figure 6 is virtually exact.

Our two different theoretical estimations for the intralayer correlation background, $G_{2\mathrm{Db}}^{++}$ and $G_{\mathrm{bg}}^{++}$, give similar results for $S(q,q_z)$ at $q_z\gtrsim 0$, which measures the transverse 2D correlation integrated over the thickness of the monolayer. For higher $q_z$, when we have to take into account the correlation structure in terms of the relative distance $z_1-z_2$, the simplest representation of the correlation background with the 2D compressibility (17(17) $\begin{array}{rl}{G}_{2\mathrm{Db}}^{++}(z_1,z_2,q)& =\int \mathrm{d}{\xi }_m\mathcal{P}({\xi }_m)G_{2\mathrm{DI}}(z_1-\xi ,z_2-\xi ,q)\\ & =h_{2\mathrm{D}}(q)\int \mathrm{d}{\xi }_m\mathcal{P}({\xi }_m)\\ & \times {\rho }_I(z_1-{\xi }_m){\rho }_I(z_2-{\xi }_m)\end{array}$(17) ) factorises that dependence as the product of the intrinsic density profiles, and the results (added to the dBW contribution), $S_{\mathrm{dBW}}^{++}(q,q_z)+S_{2\mathrm{Db}}^{++}(q,q_z)$, produce quite accurate representations of the MD data, even at large $q_z$. The optimal results are obtained when the fluctuations included in the dBW mesoscopic description are restricted to $q\le q_u\approx 2$ nm${}^{-1}$ , but this upper cutoff may be raised to $\approx 3$ nm${}^{-1}$ with only minor changes. In contrast, for the largest wavevector on the XY plane, the use of the ISM correlation background (15(15) $\begin{array}{rl}{G}_{\mathrm{bg}}^{++}(z_1,z_2,q)& \approx \int \mathrm{d}{\xi }^{\mathrm{U}}{\mathcal{P}}_{\mathrm{U}}({\xi }^{\mathrm{U}})\int \mathrm{d}{\xi }^{\mathrm{P}}{\mathcal{P}}_{\mathrm{P}}({\xi }^{\mathrm{P}})\\ & \times G_{\mathrm{Im}}(z_1-{\xi }_1^{\pm },z_2-{\xi }_2^{\pm },q),\end{array}$(15) ) $S_{\mathrm{dBW}}^{++}(q,q_z)+S_{\mathrm{bg}}^{++}(q,q_z)$, clearly underestimates $S(q,q_z)$ for $q_z\gtrsim 5$ nm${}^{-1}$ , no matter the choice of $q_u$. As SM we present similar analysis, leading to the same interpretation, for other lipid membranes and the CG MARTINI.

As for the oscillatory contribution from the interlayer correlation, the estimations of the correlation background do not add any significant contribution to the dBW mesoscopic prediction shown in Figure 6. Therefore, the strong underestimation for the amplitude of the oscillations at q = 1.34 nm${}^{-1}$, with respect to the MD results for CHARMM36-DPPC, remains unsolved.

6. Characterisation of the membrane elastic properties from the DCF

So far, we have analysed the MD results for the phosphorus atoms in the head group of the lipid molecules using two approaches: In terms of their DCF (as $G(z_1,z_2,q)$ and its Fourier transform, the structure factor $S(q,q_z)$) and in terms of their mesoscopic analysis (ISM-MD) (as the two smoothed fluctuating surfaces $z={\xi }^{\pm }(\overrightarrow{x})$). It is from that ISM analysis that we characterise the elasticity of the membrane through the functions ${\gamma }^{\mathrm{U}}(q)$, ${\gamma }^{\mathrm{P}}(q)$, ${\gamma }^{\mathrm{CU}}(q)$ and ${\gamma }^{\mathrm{m}}(q)$, directly related to the mean square fluctuations of the normal undulatory (U) and peristaltic (P) modes, or their combination as the interlayer (CU mode) and the intralayer (monolayer) fluctuations, that allow us to evaluate the mesoscopic contribution $G_{\mathrm{dBW}}(z_1,z_2,q)$ to the DCF.

In our previous work [10] we proposed and applied a deconstruction procedure for the interlayer DCF to get ${\gamma }^{\mathrm{CU}}(q)\approx \kappa q^4+\cdots$, from the density profiles ${\rho }^{\pm }(z)$ and $G^{+-}(z_1,z_2,q)$, as a method to obtain the bending modulus κ directly from the molecular positions, without having to define any mesoscopic surface. The green circles in Figure 9 show the ISM-MD result for ${\gamma }^{\mathrm{CU}}(q)$ and the green line its estimation from the direct deconstruction of the interlayer MD $G^{+-}(z_1,z_2,q)$, as obtained in our previous work [10], for the CHARMM36 representation of the DPPC tensionless bilayer membrane. The rapid increase of ${\gamma }^{\mathrm{CU}}(q)$ is produced by the decay of the correlation $⟨{\hat{\xi }}_q^{+}{\hat{\xi }}_{-q}^{-}⟩$ for $q\gtrsim 0.4$ nm${}^{-1}$ , and it was used as a practical filter for the internal fluctuations of the bilayer. Now we extend the analysis to higher q values and also to the complementary view of the membrane elasticity given by the monolayer function ${\gamma }^{\mathrm{m}}(q)$, that we estimate from the intralayer DCF.

Figure 9. The wavevector dependent surface tensions for the free DPPC All-Atoms. The circles are the values obtained from fluctuations of the intrinsic surface (ISM) shown in Figure 1, ${\gamma }^{\mathrm{m}}(q)$ (black) and ${\gamma }^{\mathrm{CU}}(q)$ (green). The lines are the ${\gamma }^{\mathrm{X}}(q)$ obtained from the deconstruction of the simulation $G(z_1,z_2,q)$. The green line ${\gamma }^{\mathrm{CU}}(q)$ of the off-diagonal quadrants $G^{+-}(z_1,z_2,q)$ obtained in our previous work [10]. The red and blue lines ${\gamma }^{\mathrm{m}}(q)$ from the deconstruction of the diagonal quadrants of bare $\mathrm{\Delta }{G}^{++}=G^{++}-G_{\mathrm{bg}}^{++}$. The red lines show the ${\gamma }^{\mathrm{m}}(q)$ obtained assuming that the correlation background is due to the compressibility, $G_{2\mathrm{Db}}^{++}$, and the blue lines are obtained using as correlation background the predicted by dBW, $G_{\mathrm{bg}}^{++}$. Dotted line with $q_{\mathrm{u}}=1.3$ nm${}^{-1}$, full line with $q_{\mathrm{u}}=1.87$ nm${}^{-1}$, and dashed line with $q_{\mathrm{u}}=2.95$ nm${}^{-1}$.

As for the analysis of the interlayer DCF in the CHARMM36 DPPC bilayer membrane, we observe in Figure 2 that the excellent agreement between the MD- ISM and the dBW deconstruction $G^{+-}(z_1,z_2,q)$ for $⟨{\hat{\xi }}_q^{+}{\hat{\xi }}_{-q}^{-}⟩\propto 1/(q^2{\gamma }^{\mathrm{CU}}(q))$ may be extended to $q\gtrsim 0.7$ nm${}^{-1}$ , when the correlation between the two lipid layers becomes negative and ${\gamma }^{\mathrm{CU}}(q)$ has a pole (not shown in Figure 9). That agreement is surprising, after the failure of the dBW structure factor in Figure 6 to reproduce the oscillations of the MD results at such large q, and it may reflect that the projection of the MD $G^{+-}(z_1,z_2,q)$ on the derivatives of the density profile, used by our deconstruction of the dBW theory, happens to filter out most of the correlation effects missed by $G_{\mathrm{dBW}}^{+-}(z_1,z_2,q)$ that we associate to the high flexibility of the lipid tails in CHARMM36, in contrast with the more rigid structures in the coarser MARTINI representation. Further studies, for other lipid molecules and other forcefield representations would be required to address this question. Another difference between the two forcefields appears when we try to represent the elasticity of the lipid membrane as $q^2/{\gamma }^{\mathrm{CU}}(q)=1/\kappa +q^2/{\kappa }_{\theta }+\mathcal{O}(q^4)$ in terms of the usual bending modulus κ and a tilt modulus ${\kappa }_{\theta }$ [10]. With MARTINI forcefield for DPPC we had shown that both κ and ${\kappa }_{\theta }$ may be obtained from (either from the fit to the MD-ISM or to the dBW DCF deconstruction) ${\gamma }^{\mathrm{CU}}(q)$. However, obtaining the tilt modulus ${\kappa }_{\theta }$ becomes complicated with CHARMM36 [10], as the strong folding of the lipid chains destroys tilt order at low q, making the tilt contribution almost insignificant by drastically reducing the q range where tilt plays a relevant role.

For the intralayer DCF, our simple proposal (16(16) $G_{2\mathrm{DI}}(z_1,z_2,q)=h_{2\mathrm{D}}(q)\sum _{\nu =\pm }{\rho }_I^{\nu }(z_1){\rho }_I^{\nu }(z_2).$(16) ) to estimate its correlation background $G_{2\mathrm{Db}}^{++}(z_1,z_2,q)$ as a 2D total correlation function $h_{2\mathrm{D}}(q)$ (with the extension in $z_{1,2}$ given by the Gaussian convolution of the intrinsic density profile) provides the opportunity to assume that the MD results for $G^{++}(z_1,z_2,q)-G_{2\mathrm{Db}}^{++}(q,z_1,z-2)$ may be taken as a $G_{\mathrm{dBW}}^{++}(z_1,z_2,q)$ from which the (unknown) function ${\gamma }^{\mathrm{m}}(q)$ may be inferred from MD results that do not require the definition and computation of any mesoscopic surface ${\xi }^{\pm }(\overrightarrow{x})$. The black circles in Figure 9 give the MD-ISM results for ${\gamma }^{\mathrm{m}}(q)$, with a non-monotonous shape reflecting the strong impact of those internal modes to describe the fluctuations of each of the two monolayers. Applying the same deconstruction method but now based on the assumption $G^{++}(z_1,z_2,q)-G_{2\mathrm{Db}}^{++}(q,z_1,z-2)\approx G_{\mathrm{dBW}}(z_1,z_2,q)$, we may reproduce quite accurately that complex shape for ${\gamma }^{\mathrm{m}}(q)$ just from the MD data of $\rho (z)$, $G(z_1,z_2,q)$ and $h_{2\mathrm{D}}(q)$, without any use of a mesoscopic (ISM) construction. The results are quite robust with respect to the choice of the wavevector cutoff $q_u$ used for the deconstruction of BW series (8(8) $\begin{array}{rl}& G_{\mathrm{dBW}}({\overrightarrow{x}}_{12},z_1,z_2)\equiv \sum _{\alpha ,\beta =\pm }{G}_{\mathrm{dBW}}^{\mathrm{\alpha \beta }}({\overrightarrow{x}}_{12},z_1,z_2)\\ & \equiv \sum _{\alpha ,\beta =\pm }(⟨{\varrho }_{\xi }^{\alpha }({\overrightarrow{x}}_1,z_1){\varrho }_{\xi }^{\beta }({\overrightarrow{x}}_2,z_2){⟩}_{\xi }-{\rho }^{\alpha }(z_1){\rho }^{\beta }(z_2))\\ & =\sum _{\alpha ,\beta =\pm }\int \int \mathrm{d}{\xi }_1^{\alpha }\mathrm{d}{\xi }_2^{\beta }{\mathcal{P}}^{\mathrm{\alpha \beta }}({\xi }_1^{\alpha },{\xi }_2^{\beta };x_{12})\\ & \times ({\rho }_{\mathrm{Im}}(\alpha (z_1-{\xi }_1^{\alpha })){\rho }_{\mathrm{Im}}(\beta (z_2-{\xi }_2^{\beta }))\\ & -{\rho }^{\alpha }(z_1){\rho }^{\beta }(z_2)\phantom{(\beta (z_2-{\xi }_2^{\beta }))})\end{array}$(8) ), and even with respect to the method used to estimate the correlation background, although (as in the comparison for $G(z_1,z_2,q)$) we get the optimal results with the simplest choice based on $h_{2\mathrm{D}}(q)$ and with any cutoff $q_u\le 2$ nm${}^{-1}$ . In the SM, we can see that the reconstruction method to obtain ${\gamma }^{\mathrm{m}}(q)$ also works, and it even performs better, for a MARTINI description of the membrane (see Figure 10 SM).

7. Discussion and conclusions

We compare two descriptions of lipid bilayer membranes, both using the same MD data for the positions of the phosphorus atom in the molecular head group of the lipids. The mesoscopic description uses the smoothed fluctuating surfaces $z={\xi }^{\pm }(\overrightarrow{x})$ for each lipid layer. The molecular descriptions uses the density correlation function (DCF), $G(z_1,z_2,q)$, Fourier transformed on the $\overrightarrow{x}=(x,y)$ plane, and as the structure factor $S(q,q_z)$ when Fourier transformed also in $z_1-z_2$. We use the same MD computer simulations of POPC and DPPC model membranes as in a previous work [10] to get direct access to $G(z_1,z_2,q)$, and the ISM method [14] to define and calculate ${\xi }^{\pm }(\overrightarrow{x})$ for the ‘upper’ $(+)$ and ‘lower’ $(-)$ monolayers of the bilayer membrane. We link these two (mesoscopic versus molecular) representations within the theoretical framework developed by Bedeaux and Weeks [11], that is extended as a double Bedeaux-Week (dBW) approach for two fluctuating surfaces [10], rather than the single one used in the Capillary Wave Theory (CWT) used in liquid surfaces. In the CWT framework [14, 23] we describe the mesoscopic elastic properties of the bilayer membrane through the wavevector dependent surface tension ${\gamma }^{\mathrm{U}}(q)$, defined in terms of the mean square Fourier component amplitudes for the fluctuating undulations (U) of the bilayer membrane as a whole, and ${\gamma }^{\mathrm{m}}(q)$ that gives the fluctuations of each lipid monolayer (i.e. one or the other edge of the membrane). The peristaltic (P) fluctuations in the local thickness of the bilayer, and the correlation between the fluctuations of the two layers (CU) are also represented as functions ${\gamma }^{\mathrm{P}}(q)$ and ${\gamma }^{\mathrm{CU}}(q)$ given by algebraic combinations of ${\gamma }^{\mathrm{U}}(q)$ and ${\gamma }^{\mathrm{m}}(q)$ [14]. The shape of these functions in Figure 9 reflects the complexity of the elastic properties of the lipid membrane over the range of q accessible to the MD simulations.

The general conclusion of this work is that the dBW formalism describes appropriately the density-density correlation and the structure factor of the simulations, at least up to $q\lesssim 0.6$ nm${}^{-1}$ and, complemented with the fluctuations due to the compressibility of each layer for the intralayer correlations, up to $q\lesssim 1.$ nm${}^{-1}$ . This theoretical technique allows a practical method to characterise the elastic properties of the lipid bilayer membrane through a set of functions ${\gamma }^{\mathrm{x}}(q)$ (for x = U, m, P or CU), that may be obtained directly from the DCF extracted by computer simulations, avoiding the need to define mesoscopic continuous fields. The usual mesoscopic elastic constants (bending modulus κ, tilt modulus ${\kappa }_{\theta }$,…) are the parameters used to characterise the set of ${\gamma }^{\mathrm{x}}(q)$ functions, restricted by their exact relationships (see Equation (2(2) $\begin{array}{rl}& \frac{1}{{\gamma }^{\mathrm{U}}(q)}+\frac{1}{{\gamma }^{\mathrm{P}}(q)}=\frac{1}{{\gamma }^{\mathrm{m}}(q)}\text{and}\\ & \frac{1}{{\gamma }^{\mathrm{U}}(q)}-\frac{1}{{\gamma }^{\mathrm{P}}(q)}=\frac{1}{{\gamma }^{\mathrm{CU}}(q)}.\end{array}$(2) )).

To reach this conclusion we have addressed two key questions; (1) What is the contribution of mesoscopic surface fluctuations to $G(z_1,z_2,q)$, if we use the ISM-MD data for ${\gamma }^{\mathrm{x}}(q)$ as input for the theoretical dBW representation? And what is the correlation background $G_{\mathrm{bg}-\mathrm{MD}}(z_1,z_2,q)=G(z_1,z_2,q)-G_{\mathrm{dBW}}(z_1,z_2,q)$ left out from that mesoscopic dBW prediction? (2) A more challenging question is: How (and how accurately) can we get the functions, ${\gamma }^{\mathrm{x}}(q)$ directly from the MD results for $G(z_1,z_2,q)$, without the need to define and use a mesoscopic representation? This is a question whose answer has a significant pragmatic consequence, as a way forward to compute elastic constants using a simpler method, that does not require to compute mathematical surfaces nor tilt fields to represent the instantaneous molecular configurations sampled along the MD simulation.

For very low wavevectors, $q\lesssim 0.25$ nm${}^{-1}$ , the large contribution from the membrane undulations to $G(z_1,z_2,q)$ dominates, and ${\gamma }^{\mathrm{U}}(q)\approx {\gamma }^{\mathrm{m}}(q)\approx {\gamma }^{\mathrm{CU}}(q)\approx {\gamma }_o+\kappa q^2+\cdots$ are given in terms of the Helfrich surface Hamiltonian [1], so that the bending modulus κ of the membrane could be directly extracted. However, for typical q values accessible in state-of-the-art MD simulations, addressing the questions above, requires the distinction of the intralayer and interlayer contributions to $G(z_1,z_2,q)$, i.e. the DCF between lipids in the same monolayer and in opposite monolayers. In a previous work [10] we considered the later, and demonstrated that it gives direct access to ${\gamma }^{\mathrm{CU}}(q)$ and to κ, since it is only weakly affected by the internal fluctuations of the bilayer, at least up to $q\lesssim 0.6$ nm${}^{-1}$.

In the present work we have addresses mainly the intralayer component of the DCF $G^{++}(z_1,z_2,q)$, associated to the function ${\gamma }^{\mathrm{m}}(q)$ which may feature a rather complex shape, particularly if we use CHARMM36 representation, as presented in Figures 1 and 9. The strong divergence between ${\gamma }^{\mathrm{m}}(q)$, ${\gamma }^{\mathrm{U}}(q)$ and ${\gamma }^{\mathrm{CU}}(q)$ for $q\gtrsim 0.4$ nm${}^{-1}$ , i.e. for wavelength $2\pi /q$ comparable to the mean membrane thickness $d\approx 5$ nm, indicates that these functions represent quite different elastic properties, which go beyond their common bending modulus. The interest of investigating the monolayer ${\gamma }^{\mathrm{m}}(q)$ comes from its interpretation as the elastic properties of the membrane edge, which would be sampled by any nanometric object (e.g. a protein approaching the membrane from the aqueous solution). Correlating, through the DCF obtained in MD simulations, the molecular characteristics of the lipid molecules with the softness or rigidity of the membrane edge, may be a very useful tool in molecular biophysics.

Moreover, we identify in our analysis the correlation background corresponding to fluctuations in the 2D density distribution (of lipid molecules head groups) on the monolayer surface, ${\xi }^{+}(\overrightarrow{x})$, rather than to fluctuations of the surface shape. For $q\gtrsim 1$ nm${}^{-1}$ the MD result for $G^{++}(z_1,z_2,q)$ cannot be recovered without that 2D correlation background contribution. The relevance (and the complexity) of the in-surface structure of biological lipid membranes is well known; but even in simple model membranes, like those studied here, their compressibility as a dense 2D liquid become dominant in the structure factor $S(q,q_z)$ at low $q_z$, since the role of the shape fluctuations vanished at $q_z=0$ (i.e. for the total integral of $G(z_1,z_2,q)$ over $z_1$ and $z_2$). Therefore, that contribution to the DCF may be most relevant to describe processes as the (partial) insertion of a membrane protein across the lipid layer. Our representation in terms of the local thickness and the local 2D density of each lipid layer contains (in an implicit way) the fluctuations that in other methods are described as a local tilt for the orientation of the lipid tails [6], since the membrane interior remains filled by the hydrophobic tails with rather low fluctuations in their 3D local density. Therefore, we may extract from $q^2/{\gamma }^{\mathrm{CU}}(q)=1/\kappa +q^2/{\kappa }_{\theta }+\mathcal{O}(q^4)$ the elastic constant ${\kappa }_{\theta }$ that in other methods is extracted from the fluctuations in the molecular tilt field.

With respect to the interlayer DFC $G^{+-}(z_1,z_2,q)$, the present work demonstrates the qualitative difference between the two different forcefield representations that we have used to model the DPPC tensionless lipid membranes. In the coarse-grained MARTINI forcefield representation, the correlation between the two lipid layers $⟨{\hat{\xi }}_q^{+}{\hat{\xi }}_{-q}^{-}⟩\equiv \mathrm{kT}/(q^2{\gamma }^{\mathrm{CU}}(q)A_o)$ decays monotonically towards zero for large q, as shown in Figure 2. In contrast, the all-atom CHARMM36 predicts that $⟨{\hat{\xi }}_q^{+}{\hat{\xi }}_{-q}^{-}⟩$ changes sign at $q\approx 0.8$ nm${}^{-1}$ , and has a minimum before raising to zero at larger q. This is reflected in Figure 1 in having ${\gamma }^{\mathrm{P}}(q)<{\gamma }^{\mathrm{U}}(q)$. These results demonstrate that for large q vectors, the CHARMM36-DPPC lipid membrane is softer with respect to modulations of the membrane thickness than membrane undulations. In contrast, for the MARTINI predicts the undulation mode is the softest one for all q (see Figure 5 of the SM ). The different fluctuations generated by all-atom (CHARMM36) or CG (MARTINI) forcefields result in a non-monotonic shape of ${\gamma }^{\mathrm{m}}(q)$, or a positive derivative of ${\gamma }^{\mathrm{m}}(q)$ for all q, respectively. Note that these differences do not affect the fit of $q^2/{\gamma }^{\mathrm{CU}}(q)=1/\kappa +q^2/{\kappa }_{\theta }+\mathcal{O}(q^4)$ at very low q to obtain the bending modulus κ [10]. However, for CHARMM36, obtaining the tilt modulus ${\kappa }_{\theta }$ from this fit becomes more problematic than for the MARTINI [10], because the folding of the lipid chains reduce the q range where tilt plays an relevant role, or even it makes uncertain the computation of a mesoscopic tilt field from the molecular configurations.

Qualitative differences between the two forcefield representation for the same lipid membrane are also evident in the MD DCFs and their mesoscopic dBW predictions. For the interlayer components with MARTINI we get $G^{+-}(z_1,z_2,q)\approx G_{\mathrm{dBW}}^{+-}(z_1,z_2,q)$ for all q. This observation supported our proposal for the deconstruction method introduced in ref. [10] to extract ${\gamma }^{\mathrm{CU}}(q)=\kappa q^2+\cdots$ from the MD result for $G^{+-}(z_1,z_2,q)$. With CHARMM36 the mesoscopic $G_{\mathrm{dBW}}^{+-}(z_1,z_2,q)$ fails to reproduce the MD results when $⟨{\hat{\xi }}_q^{+}{\hat{\xi }}_{-q}^{-}⟩$ becomes negative, for $q>0.8n{m}^{-1}$. Although our proposed method to extract the function ${\gamma }^{\mathrm{CU}}(q)$ overcomes this problem applying the procedure only within the range q>0.8  nm${}^{-1}$ , it is clear that the density correlations of all-atom (CHARMM36) bilayers do not fulfill a key CWT assumption, namely rigid local shifts of the intrinsic density profile, associated to the undulations of a fluctuating (mathematical) surface. It is logical to assume that the correlations reported at high q, using the all-atom forcefield, highlight some qualitative limitations of the coarser description with MARTINI forcefield. For instance, the different levels of detail of CG vs AA models representing the lipid bilayer are reflected in different stress distributions across lipid membranes. While early analyses of the stress profile across the bilayer show qualitatively similar profiles [12], it is known that all-atom forcefields predict more oscillations in the profile, particularly around the head group region, revealing differences in the intra/intermolecular interactions at short length scales, of the order of the CG-bead size. These differences might influence the behaviour observed in all-atom vs CG system when considering tension vs tensionless states However, it cannot be ruled out that the atomistic detail of the CHARMM36 representation introduces a peculiar behaviour for these model lipids. Resolving this problem goes beyond the main objective of this work, but we raise a word of caution about the interpretation of MD results for lipid bilayer membranes, particularly in the high q vector regime, since the results might be strongly dependent on the level of coarsening of the bilayer. Indeed, we demonstrated on our previous work [10] that the bending modulus of all-atom CHARMM36 and MARTINI forcefields are similar, but the tilt moduli, which quantifies the tilt of the lipid chains relative of the bilayer plane, are very different for all-atom or CG representations. Further work is needed to investigate the generality of these observations to other lipid bilayers.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the Spanish Secretariat for Research, Development and Innovation under: Grant No. PID2022-139776NB-C66, Grant No. PID2022-117080RB and from the Maria de Maeztu Programme for Units of Excellence in R&D [grant number CEX2018-000805-M], and by The Leverhulme Trust under UK grant RPG-2018-384. FB would like to acknowledge the Imperial College Research Computing Service, DOI: 10.14469/hpc/2232. The authors report there no completing interest to declare.