# Tight cohesion between glycolipid membranes results from balanced water-headgroup interactions

Matej Kanduc, Alexander Schlaich, Alex H. de Vries, Juliette Jouhet, Eric

Marechal, Bruno Deme, Roland R. Netz, Emanuel Schneck

► To cite this version:

Matej Kanduc, Alexander Schlaich, Alex H. de Vries, Juliette Jouhet, Eric Marechal, et al. Tight cohesion between glycolipid membranes results from balanced water-headgroup interactions. Nature Communications, 2017, 8, pp.14899. 10.1038/ncomms14899. hal-01691818

HAL Id: hal-01691818

https://hal.science/hal-01691818

Submitted on 24 Jan 2018

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

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

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Tight cohesion between glycolipid membranes results from balanced water-headgroup interactions

Matej Kanduč<fnref n="1" />,², Alexander Schlaich², Alex H. de Vries³, Juliette Jouhet⁴, Eric Maréchal⁴, Bruno Demé⁵, Roland R. Netz² & Emanuel Schneck⁶

Membrane systems that naturally occur as densely packed membrane stacks contain high amounts of glycolipids whose saccharide headgroups display multiple small electric dipoles in the form of hydroxyl groups. Experimentally, the hydration repulsion between glycolipid membranes is of much shorter range than that between zwitterionic phospholipids whose headgroups are dominated by a single large dipole. Using solvent-explicit molecular dynamics simulations, here we reproduce the experimentally observed, different pressure-versus-distance curves of phospholipid and glycolipid membrane stacks and show that the water uptake into the latter is solely driven by the hydrogen bond balance involved in non-ideal water/sugar mixing. Water structuring effects and lipid configurational perturbations, responsible for the longer-range repulsion between phospholipid membranes, are inoperative for the glycolipids. Our results explain the tight cohesion between glycolipid membranes at their swelling limit, which we here determine by neutron diffraction, and their unique interaction characteristics, which are essential for the biogenesis of photosynthetic membranes.

<footnote n="1">¹Soft Matter and Functional Materials, Helmholtz-Zentrum Berlin für Materialien und Energie, Hahn-Meitner-Platz 1, D-14109 Berlin, Germany. ²Department of Physics, Freie Universität Berlin, Arnimallee 14, D-14195 Berlin, Germany. ³Groningen Biomolecular Sciences and Biotechnology (GBB) Institute and Zernike Institute for Advanced Materials, University of Groningen, Nijenborgh 7, 9747 AG Groningen, The Netherlands. ⁴Laboratoire de Physiologie Cellulaire et Végétale, CNRS, CEA, INRA, Université Grenoble Alpes, CEA Grenoble, 17 rue des Martyrs, F-38000 Grenoble, France. ⁵Institut Laue-Langevin, 71 avenue des Martyrs, F-38042 Grenoble Cedex 9, France. ⁶Max Planck Institute of Colloids and Interfaces, Am Mühlenberg 1, D-14476 Potsdam, Germany. Correspondence and requests for materials should be addressed to M.K. (email: matej.kanduc@helmholtz-berlin.de) or to E.S. (email: schneck@mpikg.mpg.de).</footnote>
Amphiphilic lipids are the fundamental building blocks of biological membrane bilayers. Regarding the chemical structure of their hydrophilic headgroup, neutral membrane lipids can be divided into two main classes. The first are lipids with a headgroup chemistry dominated by one large electric dipole (see Fig. 1a), such as the most abundant phospholipid species phosphatidylcholine (PC). The second class involves lipids whose headgroups comprise multiple small electric dipoles, typically polar hydroxyl (OH) groups, such as glycolipids (see Fig. 1b). In nature, membranes in different cell compartments exhibit largely different lipid compositions. Highly dynamic and loosely packed membrane systems, for instance the endoplasmic reticulum or Golgi membranes, which belong to a network of endomembrane compartments all connected via vesicle budding and fusion, are rich in PC lipids¹. In contrast, structurally more steady and densely packed multilamellar membrane systems, such as myelin sheaths in vertebrates² and the photosynthetic membranes (or thylakoids) in plants³, exhibit high contents in glycolipids displaying multiple OH groups. This correlation suggests an important role of the fundamentally different headgroup architectures illustrated in Fig. 1a,b for the structural and dynamic characteristics of biological membrane systems.

Two uncharged glycolipids, mono- and di-galactosyldiacylglycerol, MGDG and DGDG, respectively, represent more than 80 per cent of the lipids in thylakoid lipid extracts⁴. The total surface area of thylakoid membranes is amplified by hierarchical organization, reaching in certain plants a total area of hundreds of square metres of thylakoids per square metre of leaves³. It is particularly striking that MGDG and DGDG are conserved from photosynthetic cyanobacteria to all chloroplasts in eukaryotes, although they are generated by completely different enzymes⁵. To what extent these glycolipids contribute to membrane stack formation and stabilization is under debate⁶. Several studies have shown that the lamellar periodicity of mature thylakoids is governed by membrane proteins⁷,⁸. However, in both cyanobacterial and eukaryotic thylakoids one finds regions in which the adjacent, glycolipid-rich membranes are in direct close contact and do not accommodate any large proteins⁹,¹⁰. Experiments on synthetic glycolipids and natural lipid extracts indicate a significant role of glycolipids in thylakoid membrane stacking: Lipid extracts from spinach chloroplasts, for instance, spontaneously form multilayers⁸. Moreover, DGDG vesicles aggregate, whereas no aggregation is observed for PC lipid vesicles¹¹. Consistently, surface force apparatus measurements revealed striking differences in the interaction between pure phospholipids and pure glycolipids (MGDG and DGDG)¹²,¹³. In a recent study on membrane stacks reconstituted from natural thylakoid lipid extracts¹⁴, it was found that water uptake significantly depends on the lipid composition. Pure DGDG membrane multilayers exhibit the strongest cohesive behaviour, that is, the weakest tendency to take up water from humid air. DGDG but also the other glycolipid mixtures were shown to swell much less than membranes composed of PC lipids¹⁵, which are enriched in the membranes of structurally more dynamic organelles. Even in excess water, the water layer between DGDG membranes remains as thin as $D_w^0 \approx 1.2$ nm, as measured in the present work by neutron diffraction (see Methods section), much less than the 2.8–3.3 nm reported for fluid PC lipid membranes¹⁶. Thus, while membrane proteins play the key functional role in thylakoids, the tight cohesion of thylakoid lipid extract membranes suggests that the lipids by themselves contribute to the tight stacking.

The interaction between membrane surfaces in water was investigated extensively in experimental and theoretical studies¹⁷,¹⁸. It is typically described on a continuum theoretical level by a superposition of electrostatic, dispersion and undulation

<figure><img src="image_3.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_6.png" /></figure>

**Figure 1 | Lipid structures and simulation set-up.** Chemical structures of a PC lipid (a) and of the glycolipid DGDG (b) as representatives of two fundamentally different lipid classes found in nature: Lipids with a headgroup chemistry dominated by one large electric dipole and lipids whose headgroups comprise multiple small electric dipoles in the form of OH groups. Both classes are schematically illustrated below the chemical structures. Dipoles are indicated by arrows. (c,d) Simulation snapshots of interacting DLPC and DGDG membranes, respectively, both at a large separation of $D_w = 2.3$ nm. With periodic boundary conditions in all three directions, the simulations represent a periodic stack of membranes with adjustable hydration level. The simulation boxes are indicated with bright rectangles. For illustration, water molecules are only shown in the lower half of the box.
forces, as well as empirical expressions for the hydration force¹⁹,²⁰. In one example, Ricoul et al.²¹ studied the interaction of mixed bilayers of cationic surfactants and glycolipids in aqueous environments. The reduced swelling of glycolipids was explained by an adhesion energy between bilayers varying linearly with the glycolipid molar fraction. However, this level of description does not account for chemical details of the interacting surfaces, which is of great importance at angstrom to nanometre separations. It is therefore now accepted that the molecular structure of surfaces and intervening solvent has to be taken into account explicitly to correctly treat the interaction and to reveal its physical mechanisms on a quantitative level on these length scales¹⁸,²². Atomistic molecular dynamics (MD) simulations allow the description of biomolecular systems at full chemical detail. However, when it comes to the interaction of extended surfaces such as membranes across an aqueous medium, the explicit treatment of water molecules has made it difficult to work at the correct water chemical potential, the latter being the key control parameter in situations with variable hydration²². Over the last few years, we have established methodology to determine the chemical potential of water and the ensuing interaction pressures in atomistic MD simulations with high precision²³,²⁴. This enables us to investigate membrane interactions on a chemically detailed and mechanistic level.

In the present work, we use atomistic MD simulations to compare the interaction mechanisms of glycolipid and phospholipid membranes. Our simulations quantitatively reproduce the experimentally observed, different pressure-versus-distance curves. Further analysis identifies the hydrogen bond balance as the driving force for the water uptake into the DGDG membranes. The associated repulsion is of short range. Water structuring effects and lipid configurational perturbation, the more long-ranged repulsion mechanisms acting between PC lipid membranes, are found to be irrelevant for DGDG.

## Results

**Area per lipid and pressure-distance curves.** The computer models of the hydrated lipid membranes employ atomistic representations of lipids and water molecules (see Fig. 1c,d for snapshots of PC lipid and DGDG membrane simulations, respectively). With periodic boundary conditions in out-of-plane and in-plane directions, the models represent infinitely extended periodic stacks of membranes in the fluid $L_\alpha$ phase with adjustable hydration level. The details of the model and simulation procedures are described in the Methods section.

Figure 2a shows the average projected area per lipid, $A_I$, of the DGDG membranes as a function of their separation $D_w$ (see Methods section for the definition of $D_w$). In the plot, open squares indicate experimental data by Shipley et al.²⁵ obtained at $T=293$ K, filled squares indicate our simulation results with semi-isotropic pressure coupling at $T=300$ K. The obtained quantitative agreement between simulation and experiment is striking. According to the low area thermal expansion coefficient of the DGDG membranes, $\alpha_A = (1.1 \pm 0.3) \times 10^{-3} \text{ K}^{-1}$ as deduced from the simulations (see Supplementary Fig. 1 and Supplementary Note 1), the temperature difference between the experiments by Shipley et al.²⁵ and our simulations affects the area per lipid by only ≈ 0.006 nm², so that it can be safely neglected in this comparison. It is seen that no significant change in $A_I$ occurs as the hydration level is varied. The simulation data are only slightly scattered around the hydration-averaged value $A_I^0 = 0.78 \text{ nm}^2$, indicated in the plot with a horizontal dashed line. This value was used for the production runs with fixed $A_I$. As shown in Supplementary Fig. 2 and described in Supplementary Note 2, variations in the force field parameters qualitatively do

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 2 | Comparison between experiments and simulations.** (a) Area per DGDG and PC lipid headgroup as functions of the water layer thickness $D_w$. Filled symbols indicate simulation results obtained with DGDG and with the PC lipid DLPC with semi-isotropic pressure coupling. Open symbols indicate experimental results for DGDG²⁵ and DLPC¹⁶ membranes. Error bars for the simulation data represent 1 s.d. of uncertainty and were estimated from the scatter of the points around the plateau value at high hydration. For DLPC, they are smaller than the size of the symbols. The experimental error can be estimated in the same way as $\le 0.02 \text{ nm}^2$. (b) Pressure-distance curves of DGDG¹⁴ and PC lipid membranes¹⁶,²⁸ as obtained in experiments (open symbols) and in the present simulations (filled symbols). Straight dashed lines in the semi-logarithmic plots indicate the best-matching exponential fits to the experimental data points for DGDG¹⁴ and to the combined experimental data sets for egg PC²⁸ and DLPC¹⁶. Error bars represent 1 s.d. of uncertainty.</figcaption></figure>

not alter this conclusion. Using area fluctuation analysis in the simulations (see Methods section), the area compressibility modulus of highly hydrated DGDG membranes was determined as $K_A = 0.35 \pm 0.10 \text{ J m}^{-2}$, slightly higher than the values typically reported for PC lipid membranes $K_A = 0.25 \text{ J m}^{-2}$ (ref. 26). In contrast, for PC lipids $A_I$ decreases by ≈ 20% on dehydration to the lowest hydration level investigated (two water molecules per lipid), as shown in Fig. 2a. Open circles indicate experimental data by Lis et al.¹⁶ for dilauroyl-PC (DLPC) obtained at $T=298$ K, filled circles again indicate our simulation results obtained at $T=300$ K (ref. 27). The good quantitative agreement between experiments and simulations for both DGDG and PC lipid membranes is a first indication of the quality of the employed computer model. The difference in area response of DGDG and PC lipids to the water content already points to distinct hydration mechanisms, as will
be corroborated below. The extent to which water molecules mediate the interaction between saccharide headgroups in the same membrane surface is investigated at high hydration ($n_w$ = 30). Each sugar headgroup on average is involved in 2.1 ± 0.1 direct hydrogen bonds (HBs, see further text for the definition) with other headgroups in the same surface and in 1.5 ± 0.1 'indirect' HBs, that is, HBs with water molecules that are simultaneously involved in a HB with other headgroups in the same surface. This result indicates that the interactions between the saccharide headgroups are water mediated to a considerable extent.

The swelling of membrane multilayers is commonly described in the form of pressure-distance curves, that is, plots of $Π$ versus $D_w$, where $Π$ is the so-called equivalent interaction pressure$^{15,28}$, which is related to the chemical potential $μ$ of water:

$$ \Pi = - \frac{\mu - \mu_0}{\nu_w^0}. \qquad (1) $$

Here, $μ_0$ and $ν_w^0$ denote the chemical potential and the partial molecular volume, respectively, of pure water in bulk. Note that $Π$ is not directly accessed in the experiment. Instead, the actual experimental control parameter is $μ$ (see Methods section). Figure 2b shows pressure-distance curves of DGDG (squares), together with those determined for PC lipids (circles) for comparison. Open symbols indicate experimental data by Demé et al.$^{14}$ ($T$ = 298 K), Lis et al.$^{16}$ ($T$ = 298 K) and Parsegian et al.$^{28}$ (room temperature), respectively. For the DGDG data$^{14}$, the errors associated with the interaction pressure were estimated by us as explained in the Methods section. For both DGDG and PC lipids, the interaction pressures are positive, that is, repulsive. In other words, work has to be performed to reduce $D_w$. The pressures decay approximately exponentially with increasing $D_w$ and reach magnitudes of several kbars for the lowest $D_w$ studied. The striking difference between DGDG and the PC lipids is the much steeper decay of the repulsive pressure in case of DGDG. An exponential fit to the experimental pressure-distance data yields decay lengths of $λ$ = 0.12 ± 0.01 nm and $λ$ = 0.27 ± 0.01 nm for DGDG and PC lipids, respectively. This difference manifests itself in the significantly different swelling limit of DGDG and PC lipids in the absence of dehydrating osmotic pressures, $D_w^0$. While $D_w^0$ is as large as ≈ 3 nm for PC lipids$^{16}$, corresponding to $n_w^0$ ≈ 30–35 water molecules per lipid, we find a swelling limit of only 1.2 nm for DGDG ($n_w^0$ ≈ 15) by neutron diffraction experiments, as described in the Methods section. Moreover, the clear difference in $λ$ evidences that the range of the repulsion is not controlled by the properties of water alone, as is often suggested in the literature$^{29}$.

The interaction pressures obtained in our simulations (filled symbols in Fig. 2b) via determination of $μ$ (see Methods section) are seen to be in remarkable quantitative agreement with the experimental data and fully reproduce the difference in the decay length. We remark that, for direct comparison with the experimental data based on equation (1), we translate $μ$ into $Π$ without accounting for the hydration dependence of the partial molecular volume of water (see below). For the water model employed in the simulations, $ν_w^0$ = 0.030 nm$^3$ for bulk water at 1 bar and 300 K. The good agreement of our simulations with the experimental data in terms of pressure-distance curves in Fig. 2b demonstrates that the mechanisms through which the hydrated adjacent membranes interact are well captured by the force fields employed in our simulations. In the following, we will analyse the simulation trajectories in detail to rationalize the characteristics of the interaction of DGDG membranes. To this end, we will highlight the differences with the interaction of PC lipid membranes.

**Interaction mechanisms.** The interaction of lipid membranes across a water layer involves a complex interplay of competing molecular interactions that collectively produce a relatively weak net repulsion. In the following, we analyse our simulations such as to identify the dominant repulsion mechanisms, keeping in mind that alternative ways to disentangle the different molecular contributions clearly exist. Figure 3 shows density profiles of water, headgroups and hydrocarbon chains perpendicular to the membrane plane for DGDG (panel a) and PC lipids (panel b)

<figure><img src="image_9.png" /><figcaption>**Figure 3 | Incorporation of the hydration water.** Density profiles of water, headgroups and hydrocarbon chains in hydrated DGDG (a) and PC lipid (b) membranes at $D_w$ = 0.23 nm. (c) Partial water volume $ν_w$ as a function of $D_w$ for DGDG and PC lipids. Error bars represent 1 s.d. of uncertainty.</figcaption></figure>
at a hydration level representative of significant repulsion, $D_w = 0.23$ nm ($\Pi = 1,200$ bar for DGDG and $\Pi = 2,100$ bar for PC lipids, see Fig. 2b). It can be seen that water molecules and headgroup moieties from both opposing membrane surfaces strongly overlap at the midplane, that is, in the centre of the water distribution. Water uptake will therefore affect lipid-lipid (ll), lipid-water (lw) and water-water (ww) interactions, and the interaction balance will depend on the way additional water molecules are accommodated in the layer of partially hydrated headgroups. Important insight into this can be obtained from the partial water volume, $v_w = (\partial V/\partial N_w)_{p,T}$ defined as the change of the system volume on inserting a water molecule at constant pressure and temperature. In our simulations, we evaluate $v_w$ by fitting the system volume $V(N_w)$ at various hydration levels $N_w$ and taking the derivative with respect to $N_w$. In Fig. 3c, $v_w$ is shown as a function of $D_w$ for DGDG and PC lipids. It can be seen that $v_w$ at large $D_w$ coincides with the bulk water value $v_w^0$ for both DGDG and PC lipids. While for PC lipids, this holds true for the entire hydration range, for DGDG we observe a significant deviation towards lower values, $v_w < v_w^0$, for small $D_w$. This means that on addition of a water molecule at low hydration, the system expands less than a bulk water system does. This observation reflects the non-ideal mixing of water and saccharides. The significant hydration dependence of $v_w$ for DGDG therefore shows that the use of equation (1), assuming $v_w = v_w^0$, is not always justified for the calculation of $\Pi$, neither in experiments nor in simulations. According to the relation between $\Pi$ and the free energy $G$ per area $A$ at constant ambient pressure $p$ and temperature $T$,

$$ \Pi = - \frac{1}{A} \left( \frac{dG}{dD_w} \right)_{p,T}, \quad (2) $$

the significant repulsion ($\Pi > 0$) in the short separation range indicates that the incorporation of water into the partially hydrated saccharide layer lowers the free energy. Since both water and saccharides display OH groups at high densities, the key contribution to this free energy reduction on hydration is expected to be related to the formation of HBs. On a simplified level, HBs are electrostatic interactions involving two atoms carrying negative partial charges and a hydrogen atom that is covalently bound to one of them and carries a positive partial charge. To extract HB numbers from our simulation trajectories, we employ the widely used Luzar-Chandler geometric criterion³⁰, according to which a HB is present if the distance between donor and acceptor atoms is smaller than 0.35 nm and the hydrogen-donor-acceptor angle is smaller than 30°.

The total number of HBs per lipid molecule, $n_{\text{tot}}^{\text{HB}}$, is defined as the sum of ll, lw and excess ww HBs:

$$ n_{\text{tot}}^{\text{HB}} = n_{\text{ll}}^{\text{HB}} + n_{\text{lw}}^{\text{HB}} + n_{\text{ww}}^{\text{HB,ex}}, \quad (3) $$

where $n_{\text{ww}}^{\text{HB,ex}} = n_{\text{ww}}^{\text{HB}} - n_w n_{\text{bw}}^{\text{HB}}$ and $n_{\text{bw}}^{\text{HB}} = 1.796$ is the number of HBs per water molecule in bulk water, obtained in separate simulations. Figure 4a shows the change in $n_{\text{tot}}^{\text{HB}}$ for DGDG on dehydration, $\Delta n_{\text{tot}}^{\text{HB}}(D_w) = n_{\text{tot}}^{\text{HB}}(D_w) - n_{\text{tot}}^{\text{HB}}(\infty)$. It can be seen that $\Delta n_{\text{tot}}^{\text{HB}}$ increases on swelling, that is, with rising $D_w$, until it reaches zero at large $D_w$. This means that water uptake increases the overall number of HBs. In fact, the decay lengths of $\Pi$ and $n_{\text{tot}}^{\text{HB}}$ coincide ($\lambda = 0.12$ nm for both, see Figs 2b and 4a where the data are compared with exponential fits), suggesting that hydrogen bonding is at the heart of the swelling mechanism of DGDG membranes. A more rigorous test of this hypothesis is presented in Fig. 4b, where $\Delta n_{\text{tot}}^{\text{HB}}$ is plotted versus the interaction free energy per lipid $G/N_l$ obtained by integrating the pressure-distance curve, $G = -A \int \Pi dD_w$ according to equation (2). The plot reveals a linear relation between hydrogen bond number $\Delta n_{\text{tot}}^{\text{HB}}$ and free energy $G$, with zero intercept. The linear fit yields

the slope $d(G/N_l)/d(\Delta n_{\text{tot}}^{\text{HB}}) = -21 \pm 3$ kJ mol$^{-1}$, which is comparable to experimental estimates of HB free energies³¹. In other words, the repulsion between DGDG membranes can be described entirely in terms of the HB balance. In the inset of Fig. 4a, we further show that with increasing water uptake the number of lw HBs grows at the expense of ll and ww HBs. However, the exponential decay length for all three curves, $\lambda = 0.26 \pm 0.01$ nm, is not reflected in the pressure-distance curve of DGDG (indicated by Fig. 2b), indicating that ll, lw and ww interactions balance each other and thus compensate each other almost perfectly. This is no surprise in view of the similar chemistry of water and saccharide OH groups. For PC lipids, the situation is fundamentally different: the chemistry of the PC headgroup, dominated by a large single dipole of which the positive charge is encased by three hydrophobic methyl groups (see Fig. 1a), has little in common with that of water. In fact, it has been shown that ll, lw and ww interactions for PC lipid membranes are disparate and an energetic preference for lw interactions gives rise to strong repulsion between PC lipid membranes at low hydration²⁴. As shown in Fig. 4c, there is no simple proportionality between $G$ and $\Delta n_{\text{tot}}^{\text{HB}}$ for PC lipids (see Supplementary Fig. 3 and Supplementary Note 3 for the separation dependence of $\Delta n_{\text{tot}}^{\text{HB}}$). Moreover, in the limit of high hydration (shaded in light blue in Fig. 4c) $G$ varies while the HB number has already saturated ($\Delta n_{\text{tot}}^{\text{HB}} \approx 0$). This behaviour reflects more long-ranged, non-HB-related repulsion mechanisms between PC lipid membranes at higher hydration, giving rise to the larger decay length in the pressure-distance curves in Fig. 2b.

Two mechanisms are responsible for this long-ranged repulsion: water structuring effects via the strong orientational polarization of the water layers interacting with the lipid headgroups³², and the configurational entropy of the lipids³³. In the following, we show that both mechanisms are essentially inoperative for glycolipids like DGDG. Figure 5a presents density profiles of water and lipids in hydrated DGDG and PC lipid membranes at a large separation of $D_w = 2.3$ nm. It is seen that the water density profiles are nearly identical for the two lipid types. In contrast, as shown in Fig. 5b, DGDG and PC lipid membranes lead to very different profiles of water orientation perpendicular to the membrane surfaces, $\langle \cos \theta_w \rangle$ (see inset for the definition of $\theta_w$). For the PC lipids, the water dipoles close to the membrane surfaces at $z = \pm D_w/2 = \pm 1.15$ nm (see vertical dashed lines) are strongly oriented and significant orientation extends virtually all the way to the centre of the water layer (see zoom-in in Supplementary Fig. 4 and Supplementary Note 4), where it has to vanish by symmetry. For DGDG, on the other hand, water orientation is pronounced only in the poorly hydrated inner headgroup regions and insignificant inside the water layer. This dissimilar behaviour can be understood from the different headgroup structures illustrated in Fig. 1a,b: The large and directionally correlated single electric dipoles of PC lipids induce strong water orientation, while the rather isotropically oriented OH groups of the glycolipids do not. Significant repulsion due to water structuring can thus not be expected for DGDG membranes. We move on with Fig. 5c, which shows the angular distributions of DGDG and PC lipid headgroups with respect to the membrane normal for large ($D_w = 2.3$ nm, solid lines) and small ($D_w = 0.6$ nm, dashed lines) separations. Headgroup angles $\theta_l$ are defined by vectors between pairs of headgroup atoms, as illustrated in the figure inset: P-N for PC lipids and opposing headgroup carbons $C_1-C_2$ for DGDG. The distributions are weighted with a factor $1/\sin(\theta_l)$, so that a constant distribution would correspond to a random orientation. It can be seen that the angular distribution of the PC headgroups is quite broad for large separations but narrows considerably on dehydration. In fact, more upright conformations with $\theta_l$ smaller
<figure><img src="image_10.png" /><figcaption>**Figure 4** | **Hydrogen bonds.** (a) Total number of HBs per DGDG lipid as a function of the membrane separation. The dashed line indicates an exponential fit with decay length $\lambda = 0.12$ nm. Inset: II, Iw and ww HBs. In the plot, ww refers to the excess HB number with respect to bulk water, $n_{ww}^{HB} - n_w n_{bw}^{HB}$, see main text. Dashed lines indicate exponential fits with the same decay length $\lambda = 0.26 \pm 0.01$ nm for all three curves. (b) The total number of HBs versus the interaction free energy for DGDG. Dashed line: Linear regression through all data points. (c) The same plot for PC lipids. Shaded region: High hydration limit where $\Delta n_{tot}^{HB} \approx 0$. The error bars in all three panels represent 1 s.d. of uncertainty.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 5** | **Origin of the long-range repulsion.** (a) Density profiles of water and lipids in hydrated DGDG and PC lipid membranes at a large separation, $D_w = 2.3$ nm. Vertical lines indicate the membrane surfaces at $\pm D_w/2$. (b) Water orientation profiles $\langle \cos \theta_w \rangle$ at the same separation. Inset: Definition of the water dipole angle $\theta_w$. (c) Distributions of DGDG and PC lipid headgroup orientations with respect to the membrane normal for large ($D_w = 2.3$ nm, solid lines) and small ($D_w = 0.6$ nm, dashed lines) separations. Inset: Definition of the headgroup vectors.</figcaption></figure>

than approximately 45° are largely suppressed for $D_w = 0.6$ nm. In contrast, the headgroup angular distribution for DGDG is much narrower and, importantly, virtually unaffected by dehydration. (A comparison at the same interaction pressure instead of the same separation yields the same picture, see Supplementary Fig. 5 and Supplementary Note 5.) The comparison in Fig. 5c indicates that entropic repulsion due to a lipid configurational perturbation is of minor relevance for the glycolipid membranes. The latter conclusion is further corroborated by decomposing the interaction free energy G into the enthalpic and entropic contributions, which we show in Supplementary Fig. 6 and discuss in Supplementary Note 6.

## Discussion

With the preceding analysis, we have shown that the water uptake into DGDG membranes is solely driven by the HB balance involved in non-ideal water/sugar mixing. The ensuing repulsion is of much shorter range than that between phospholipid membranes, which exhibit more long-ranged swelling mechanisms. The equilibrium separation between membranes in excess water, $D_w^0$, is known to be governed by the balance between repulsive hydration forces, as quantified in the present study, and the ubiquitous van der Waals (vdW) attraction between the membranes across the water layer<sup>19,20</sup>. More short-ranged repulsion for DGDG thus coincides with smaller $D_w^0$. Because the vdW attractive potential obeys a $D_w^{-2}$ scaling<sup>34</sup>, the depth of the adhesive free energy minimum at the equilibrium separation, $G_{adh} \equiv G(D_w^0)$ (which we approximate as the vdW minimum neglecting the tail of the exponential repulsion), scales as $(D_w^0)^{-2}$ and is thus significantly larger for the glycolipids. Comparing DGDG ($D_w^0 = 1.2$ nm, measured by neutron diffraction in the

<figure><img src="image_12.png" /><figcaption>**Figure 6** | **Dominant species of DGDG.** Chemical structures of dalDGDG (top) and palDGDG (bottom) lipids, differing in the saturation of the hydrocarbon chains.</figcaption></figure>

present work) and PC lipids ($D_w^0 \approx 3$ nm), we numerically obtain for the ratio of the adhesion strengths $G_{adh}^{DGDG}/G_{adh}^{PC} \propto (D_w^{0,DGDG}/D_w^{0,PC})^{-2} \approx (1.2 \text{ nm}/3.0 \text{ nm})^{-2} \approx 6$, when neglecting differences in the Hamaker constants. With that, our analysis provides the physical explanation for the stronger cohesion of glycolipid membranes and their tendency to form stable lamellar structures. The striking differences in the interaction mechanisms of DGDG and PC lipids, as representatives of the two fundamentally different lipid classes defined in Fig. 1a,b, underline the impact the headgroup design has on membrane interactions. In fact, for certain saccharides, strong attractive contributions to the interfacial force balance have been reported<sup>20,35</sup>. The characteristics of the repulsion between lipid membranes (or hydrophilic surfaces in general) in terms of strength, range and mechanisms are not governed by a universal, water-inherent mechanism, but highly sensitive to the headgroup chemistry.
<figure><img src="image_13.png" /><figcaption>**Figure 7** | Lamellar period at excess hydration. Neutron diffraction of pure DGDG measured in excess D₂O between two silicon wafers. (a) 2D diffraction pattern. (b) 1D reduced curve. The insert is a plot of the Bragg reflection (h=1 to 3) positions versus h. The solid line is a linear fit to the experimental points yielding a period of 54.0 ± 0.1 Å.</figcaption></figure>

Clearly, ions are expected to have significant influence on the hydration repulsion, as they affect the balance between effective water-water and water-headgroup interactions, notably HBs³⁶. In the future, it would be interesting to compare the influence of ions on the swelling of phospholipids and glycolipids.

The differences in the interaction between glycolipid and phospholipid membranes may be of relevance for protein-free regions in thylakoids, which were reported to exhibit membrane separations in a broad range of 2–4 nm (refs 37,38). In this distance range, phospholipid membranes still exhibit significant hydration repulsion, in contrast to glycolipids, for which the hydration repulsion has basically decayed to zero. This might suggest that protein-protein interactions are not the only factor influencing thylakoid architecture and we speculate that membrane interactions may therefore be relevant even for the evolution of lipid headgroup chemistry.

Complementary information on the tightly packed domains in thylakoids in terms of membrane bending rigidities may be obtained in future studies exploiting the off-specular scattering of neutrons from aligned membrane multilayers²⁰,³⁹. Finally, sugar-water interactions are not only relevant for glycolipids but also in a much broader context ranging from sugar solutions⁴⁰ and sugar surfactants⁴¹,⁴² to glycoproteins and sugar-based biomaterials⁴³. Some of the concepts presented here may also apply to those problems.

## Methods

**Set up of the computer model.** The computer model of the hydrated glycolipid bilayers employs atomistic representations of lipid and water molecules. Each of the two bilayer leaflets contains 50 DGDG molecules, see Fig. 1b. In the scattering experiments by Demé et al.¹⁴, the DGDG sample was dominated by two main DGDG populations with different acyl chains, see Fig. 6: dalDGDG with two threefold-unsaturated C18 chains (18:3) and palDGDG with one 18:3 chain at the sn-2 position and one saturated C16 chain (16:0) at the sn-1 position. In order for the simulations to be comparable to the experiments, dalDGDG and palDGDG were also realized in the computer model and mixed in the approximate ratios reported in the experiments, 80% dalDGDG and 20% palDGDG. The systems were prepared by using the insane.py script⁴⁴ available from the Martini website (www.cgmartini.nl). The script can be used to prepare coarse-grained (CG) models of bilayers of varying composition by using lipid templates. After minimization and

a short equilibration run at constant area, atomistic details were generated from the CG structures using the 'backward.py' script⁴⁵, also available from the Martini website. Mapping files from the CG to the atomistic structures were written for the lipids relevant to these systems. The atomistic structures after backmapping were energy minimized and briefly equilibrated using stochastic dynamics with very short time step⁴⁶,⁴⁷.

**Computer simulations.** All atomistic MD simulations are performed using the GROMACS simulation package⁴⁸. Simulation run files are available as Supplementary Software 1. We use the simple point charge/extended water model⁴⁹. Electrostatics is treated using the particle-mesh-Ewald method⁵⁰,⁵¹ with a 0.9 nm real-space cutoff. The Lennard-Jones (LJ) interactions are cut off at $r_{LJ} = 0.9$ nm. The simulations are performed with an integration time step of 2 fs in the canonical constant pressure ensemble with periodic boundary conditions. Temperature was maintained at T = 300 K, using the Berendsen thermostat⁵² with a time constant of 1 ps. The pressure was maintained at p = 1 bar using the Berendsen barostat with a time constant of 1 ps.

For the PC lipids, we use the united-atom Berger force field⁵³. The bilayer is composed of 72 DLPC molecules (36 in each layer) and subject to semi-isotropic pressure coupling at 1 bar (independent lateral and perpendicular simulation box scaling). For each hydration level, the membranes are equilibrated for at least 5 ns before the production runs, which have durations of 75 ns. The thermodynamic integration (TI) simulations needed for the determination of the water chemical potential (see below) has a total duration of 800 ns per hydration level.

For the glycolipid simulations, we start from the above-described initial configurations. The DGDG force field is based on GROMOS 53a6, which employs united-atom treatment of non-polar hydrogens⁵⁴ and is described in refs 46,47. Due to their pronounced hydrogen-bonding capabilities, DGDG bilayers exhibit much longer relaxation times than the PC lipid bilayers. To nevertheless achieve sufficient sampling statistics, additional measures are taken: At first, for each hydration level an initial 80 ns simulation run with semi-isotropic pressure coupling at 1 bar is performed to determine the corresponding equilibrium value of $A_i$. In the next step, however, the lateral box extension, and thus $A_i$, is then kept at the constant value $A_i^0$, as explained in the text. This treatment reduces fluctuations of the box vectors and the hydration water layer thickness and facilitates measurements of the water chemical potential. Second, for each hydration level, 15 independent sets of simulations are averaged, amounting to a total duration per hydration level of 1,000 ns of production runs and 7,500 ns of simulations for TI.

The water layer thickness (or membrane separation) $D_w$ is defined via the number of water molecules per lipid $n_w$, as $D_w = 2n_w v_w^0/A_i$, where $v_w^0 = 0.03$ nm³ is the volume of a water molecule in bulk and $A_i$ is the average area per lipid. This definition is commonly used in the experimental literature²⁸.

The equivalent interaction pressure, computed via equation (1), is in experiments obtained by controlling the water chemical potential $\mu$. In the experimental work on DGDG¹⁴, for instance, $\mu$ is adjusted via the relative humidity $h_{rel}$, as $\mu = \mu_0 + k_B T \ln(h_{rel})$.

We have recently established the methodology to determine the chemical potential of water in atomistic MD simulations with high precision. $\Pi$ is then obtained from the shift in the chemical potential in analogy with equation (1). To determine $\mu = \mu^{id} + \mu^{ex}$, we independently measure its excess and ideal contributions, $\mu^{ex}$ and $\mu^{id} = k_B T \ln(\rho)$, where $\rho$ is the water density. While $\mu$ by definition is constant over the simulation volume in thermal equilibrium, $\mu^{ex}$ and $\mu^{id}$ are not. Due to the inhomogeneous water distribution perpendicular to the membrane surface, $\mu^{ex}(z)$ and $\mu^{id}(z)$ via $\rho(z)$ are functions of the perpendicular coordinate, $z$. As a consequence, $\mu^{ex}$ and $\mu^{id}$ have to be evaluated at the same $z$ position in the simulation box to determine $\mu$. This is facilitated for minimal box vector fluctuations and thus motivates constraining $A_i$ in the production runs. We

<table><caption>**Table 1** | Bragg peak positions.</caption><thead><tr><td>Bragg order</td><td>q<sub>x</sub> (Å<sup>-1</sup>)</td><td>δq<sub>x</sub> (Å<sup>-1</sup>)</td></tr></thead><tbody><tr><td>1</td><td>0.11697</td><td>0.000130</td></tr><tr><td>2</td><td>0.23346</td><td>0.000103</td></tr><tr><td>3</td><td>0.34957</td><td>0.000302</td></tr></tbody></table>
then conveniently determine $\mu^{id}$ from $\rho$, whereas for the determination of $\mu^{ex}$ we use a computationally efficient combination of two approaches, the test particle insertion⁵⁵ and the Ti⁵⁶. Additional methodological details are described elsewhere⁵⁷. At this point, we remark that for consistent evaluations of the LJ energies, the LJ potentials should be shifted by an offset to reach zero at their cutoff distance $r_{1LJ}$. The reason is that the particle trajectories in MD are produced by computing the forces from potentials, and therefore any potential that is cut off effectively behaves as being shifted. This fact was ignored in refs 24,27, which led to slightly different pressure-distance curves.

**Area compressibility modulus.** In the simulations with unconstrained lateral box area A, the area compressibility modulus $K_A$ is calculated from the fluctuations in the area according to ref. 58

$$K_A = \frac{k_B T \langle A \rangle}{\langle (A - \langle A \rangle)^2 \rangle} \quad (4)$$

To give an error estimate for $K_A$, the total simulation time is divided into 5–10 blocks (leading to various lengths of 100–200 ns) and the value for $K_A$ is calculated for each block. The analysis was performed for hydration levels of 5 and 20 waters per DGDG lipid and yielded the following results:

$$n_w = 5; K_A = 0.85 \pm 0.20 \text{ J m}^{-2}$$

$$n_w = 20; K_A = 0.35 \pm 0.10 \text{ J m}^{-2}$$

**Experimental errors of DGDG interaction pressures.** In the original pressure-distance measurements on DGDG membranes by Demé et al.¹⁴, the errors associated with the interaction pressure were not estimated. In their work, the hydrated membrane multilayers were in equilibrium with water vapour of known chemical potential $\mu$, which varies with relative humidity $h_{rel}$ as $\mu = \mu_0 + k_B T \ln h_{rel}$. The equivalent osmotic pressure was then calculated using equation (1), viz.

$$\Pi = - \frac{k_B T}{\nu_w^0} \ln h_{rel}. \quad (5)$$

The accuracy of the measured osmotic pressure $\Pi$ therefore depends on the accuracy of the measured relative humidity $h_{rel}$. The accuracy of the freshly calibrated humidity sensor can be estimated as $\pm 2$ per cent ($\delta h_{rel} \approx 0.02$). By differentiating equation (5), we obtain the error estimate for the osmotic pressure as $\delta\Pi = (k_B T/\nu_w^0)(\delta h_{rel}/h_{rel})$. Expanding the relative humidity $h_{rel}$ in terms of the osmotic pressure $\Pi$ using equation (5) then leads to

$$\delta\Pi = \frac{k_B T}{\nu_w^0} \delta h \exp\left(\frac{\nu_w^0 \Pi}{k_B T}\right). \quad (6)$$

This expression enables us to estimate the error bars for the DGDG data in Fig. 2b.

**Neutron diffraction experiments.** The lamellar periodicity of DGDG at full hydration is measured by exposing a solid-supported, stacked DGDG sample to excess water (in the form of D₂O) between two silicon wafers (Si-Mat, Kaufering, Germany) in a sandwich configuration³⁹. Under these conditions, the lamellar stack is at maximum swelling and zero osmotic pressure. The sample has the same composition as the one used in ref. 14 in which the experimental DGDG pressure-distance curve is reported: $\approx 80\%$ dalDGDG and $\approx 20\%$ palDGDG, where dalDGDG has two threefold-unsaturated C18 chains (18:3) and palDGDG has one 18:3 chain at the sn-2 position and one saturated C16 chain (16:0) at the sn-1 position, see Fig. 6.

Neutron diffraction experiments are carried out as described in ref. 14 on the D16 instrument at the ILL Grenoble, France, using a wavelength $\lambda = 4.75$ Å and $\Delta\lambda/\lambda = 0.01$. Data analysis is performed using the ILL in-house LAMP software (www.ill.eu/instruments-support/computing-for-science/cs-software/all-software/lamp): The intensity collected on the 2D detector (Fig. 7a) is reduced to 1D by vertical integration after solid angle and detector pixel efficiency correction. The data are corrected for background by subtraction of the empty chamber signal. The lamellar periodicity (D) is determined in a linear fit to three Bragg peak positions (Fig. 7b) versus diffraction order $h$ of the lamellar ($h00$) reflections, according to Bragg's law: $q_s^{(h)} = 2\pi h/D$, where $q_s^{(h)}$ is the magnitude of the scattering vector. Each peak is fitted with a Gaussian function on top of a constant background, yielding the following $q_s$-values and errors ($\delta q_s$) for the peak position (Table 1).

The peak widths are governed by the instrument settings (notably by the collimation of the neutron beam) and do not reflect any disorder of the sample. The linear fit to the peak positions yields $D = 54.0 \pm 0.1$ Å. Assuming a bilayer thickness of 41.7 Å (ref. 25) consistently with ref. 14, this lamellar period corresponds to a bilayer separation distance of 12.3 Å.

**Data availability.** All data supporting the findings of this study are available within the article and its Supplementary Information files.

## References

1. Van Meer, G., Voelker, D. R. & Feigenson, G. W. Membrane lipids: where they are and how they behave. *Nat. Rev. Mol. Cell Biol.* **9**, 112–124 (2008).

2. Stoffel, W. & Bosio, A. Myelin glycolipids and their functions. *Curr. Opin. Neurobiol.* **7**, 654–661 (1997).

3. Boudière, L. et al. Glycerolipids in photosynthesis: composition, synthesis and trafficking. *Biochim. Biophys. Acta* **1837**, 470–480 (2014).

4. Block, M. A., Dorne, A.-J., Joyard, J. & Douce, R. Preparation and characterization of membrane fractions enriched in outer and inner envelope membranes from spinach chloroplasts. *ii. biochemical characterization*. *J. Biol. Chem.* **258**, 13287–13286 (1983).

5. Petroutsos, D. et al. Evolution of galactoglycerolipid biosynthetic pathways–from cyanobacteria to primary plastids and from primary to secondary plastids. *Prog. Lipid. Res.* **54**, 68–85 (2014).

6. Chow, W. S., Kim, E.-H., Horton, P. & Anderson, J. M. Granal stacking of thylakoid membranes in higher plant chloroplasts: the physicochemical forces at work and the functional consequences that ensue. *Photochem. Photobiol. Sci.* **4**, 1081–1090 (2005).

7. McDonnel, A. & Staehelin, L. A. Adhesion between liposomes mediated by the chlorophyll a/b light-harvesting complex isolated from chloroplast membranes. *J. Cell Biol.* **84**, 40–56 (1980).

8. Ryrif, I. J., Anderson, J. M. & Goodchild, D. J. The role of the light-harvesting chlorophyll a/b-protein complex in chloroplast membrane stacking. *Eur. J. Biochem.* **107**, 345–354 (1980).

9. Nagy, G. et al. Reversible membrane reorganizations during photosynthesis *in vivo*: revealed by small-angle neutron scattering. *Biochem. J.* **436**, 225–230 (2011).

10. Stingaciu, L.-R. et al. Revealing the dynamics of thylakoid membranes in living cyanobacterial cells. *Sci. Rep.* **6**, 19627 (2016).

11. Webb, M. S. & Green, B. R. Effects of neutral and anionic lipids on digitalactosyldiacylglycerol vesicle aggregation. *Biochim. Biophys. Acta* **1030**, 231–237 (1990).

12. Marra, J. Controlled deposition of lipid monolayers and bilayers onto mica and direct force measurements between galactolipid bilayers in aqueous solutions. *Adv. Colloid Interface Sci.* **107**, 446–458 (1985).

13. Marra, J. Direct measurements of attractive van der waals and adhesion forces between uncharged lipid bilayers in aqueous solutions. *Adv. Colloid Interface Sci.* **109**, 11–20 (1986).

14. Demé, B., Cataye, C., Block, M. A., Maréchal, E. & Jouhet, J. Contribution of galactoglycerolipids to the 3-dimensional architecture of thylakoids. *FASEB J.* **28**, 3373–3383 (2014).

15. LeNeveu, D. M., Rand, R. P. & Parsegian, V. A. Measurement of forces between lecithin bilayers. *Nature* **259**, 601–603 (1976).

16. Lis, L. J., McAlister, M., Fuller, N., Rand, R. P. & Parsegian, V. A. Interactions between neutral phospholipid bilayer membranes. *Biophys. J.* **37**, 657–665 (1982).

17. Lipowsky, R. & Sackmann, E. *Structure and Dynamics of Membranes: I. From Cells to Vesicles/II. Generic and Specific Interactions* (Elsevier, 1995).

18. Rand, R. P. & Parsegian, V. A. Hydration forces between phospholipid bilayers. *Biochim. Biophys. Acta* **988**, 351–376 (1989).

19. Leontidis, E., Aroti, A., Belloni, L., Dubois, M. & Zemb, T. Effects of monovalent anions of the hofmeister series on DPPC lipid bilayers part II: modeling the perpendicular and lateral equation-of-state. *Biophys. J.* **93**, 1591–1607 (2007).

20. Schneck, E., Demé, B., Gege, C. & Tanaka, M. Membrane adhesion via homophilic saccharide-saccharide interactions investigated by neutron scattering. *Biophys. J.* **100**, 2151–2159 (2011).

21. Ricoul, F. et al. Phase equilibria and equation of state of a mixed cationic surfactant-glycolipid lamellar system. *Langmuir* **14**, 2645–2655 (1998).

22. Schneck, E. & Netz, R. R. From simple surface models to lipid membranes: universal aspects of the hydration interaction from solvent-explicit simulations. *Curr. Opin. Colloid Interface Sci.* **16**, 607–611 (2011).

23. Kanduč, M., Schneck, E. & Netz, R. R. Attraction between hydrated phospholipic surfaces. *Chem. Phys. Lett.* **610**, 375–380 (2014).

24. Schneck, E., Sedlmeier, F. & Netz, R. R. Hydration repulsion between biomembranes results from an interplay of dehydration and depolarization. *Proc. Natl Acad. Sci. USA* **109**, 14405–14409 (2012).

25. Shipley, G. G., Green, J. P. & Nichols, B. W. The phase behavior of monogalactosyl, digitalactosyl, and sulphoquinovosyl diglycerides. *Biochim. Biophys. Acta* **311**, 531–544 (1973).

26. Rawicz, W., Olbrich, K. C., McIntosh, T., Needham, D. & Evans, E. Effect of chain length and unsaturation on elasticity of lipid bilayers. *Biophys. J.* **79**, 328–339 (2000).

27. Kanduč, M., Schneck, E. & Netz, R. R. Hydration interaction between phospholipid membranes: insight into different measurement ensembles from atomistic molecular dynamics simulations. *Langmuir* **29**, 9126–9137 (2013).

28. Parsegian, V. A., Fuller, N. & Rand, R. P. Measured work of deformation and repulsion of lecithin bilayers. *Proc. Natl Acad. Sci. USA* **76**, 2750–2754 (1979).

29. Parsegian, V. A. & Zemb, T. Hydration forces: observations, explanations, expectations, questions. *Curr. Opin. Colloid Interface Sci.* **16**, 618–624 (2011).

30. Luzar, A. & Chandler, D. Hydrogen-bond kinetics in liquid water. *Nature* **379**, 55–57 (1996).
31. Fersht, A. R. et al. Hydrogen bonding and biological specificity analyzed by protein engineering. *Nature* **314**, 235–238 (1985).

32. Marčelja, S. & Radic, N. Repulsion of interfaces due to boundary water. *Chem. Phys. Lett.* **42**, 129–130 (1976).

33. Israelachvili, J. N. & Wennerström, H. Entropic forces between amphiphilic surfaces in liquids. *J. Phys. Chem.* **96**, 520–531 (1992).

34. Parsegian, V. A. *Van der Waals Forces: A Handbook for Biologists, Chemists, Engineers, and Physicists* (Cambridge University Press, 2005).

35. Yu, Z. W., Calvert, T. L. & Leckband, D. Molecular forces between membranes displaying neutral glycosphingolipids: evidence for carbohydrate attraction. *Biochemistry* **37**, 1540–1550 (1998).

36. Marčelja, S. Hydration forces near charged interfaces in terms of effective ion potentials. *Curr. Opin. Colloid Interface Sci.* **16**, 579–583 (2011).

37. Daum, B., Nicastro, D., Austin, J., McIntosh, J. R. & Kühlbrandt, W. Arrangement of photosystem II and ATP synthase in chloroplast membranes of spinach and pea. *Plant Cell* **22**, 1299–1312 (2010).

38. Dekker, J. P. & Boekema, E. J. Supramolecular organization of thylakoid membrane proteins in green plants. *Biochim. Biophys. Acta* **1706**, 12–39 (2005).

39. Schneck, E. et al. Modulation of intermembrane interaction and bending rigidity of biomembrane models via carbohydrates investigated by specular and off-specular neutron scattering. *Phys. Rev. E* **78**, 061924 (2008).

40. Sauter, J. & Grafmüller, A. Solution properties of hemicellulose polysaccharides with four common carbohydrate force fields. *J. Chem. Theory Comput.* **11**, 1765–1774 (2015).

41. Johnsson, M. & Engberts, J. B. F. N. Novel sugar-based gemini surfactants: aggregation properties in aqueous solution. *J. Phys. Org. Chem.* **17**, 934–944 (2004).

42. Stubenrauch, C. Sugar surfactants–aggregation, interfacial, and adsorption phenomena. *Curr. Opin. Colloid Interface Sci.* **6**, 160–170 (2005).

43. Bertinetti, L., Fratzl, P. & Zemb, T. Chemical, colloidal and mechanical contributions to the state of water in wood cell walls. *New J. Phys.* **18**, 083048 (2016).

44. Wassenaar, T. A., Ingólfsson, H. I., Böckmann, R. A., Tieleman, D. P. & Marrink, S. J. Computational lipidomics with insane: a versatile tool for generating custom membranes for molecular simulations. *J. Chem. Theory Comput.* **11**, 2144–2155 (2015).

45. Wassenaar, T. A., Pluhackova, K., Böckmann, R. A., Marrink, S. J. & Tieleman, D. P. Going backward: a flexible geometric approach to reverse transformation from coarse grained to atomistic models. *J. Chem. Theory Comput.* **10**, 676–690 (2014).

46. Lopez, C. A., Sovova, Z., van Eerden, F. J., de Vries, A. H. & Marrink, S. J. Martini force field parameters for glycolipids. *J. Chem. Theory Comput.* **9**, 1694–1708 (2013).

47. van Eerden, F. J., de Jong, D. H., de Vries, A. H., Wassenaar, T. A. & Marrink, S. J. Characterization of thylakoid lipid membranes from cyanobacteria and higher plants by molecular dynamics simulations. *Biochim. Biophys. Acta* **1848**, 1319–1330 (2015).

48. Van Der Spoel, D. et al. GROMACS: fast, flexible, and free. *J. Comput. Chem.* **26**, 1701–1718 (2005).

49. Berendsen, H. J. C., Grigera, J. R. & Straatsma, T. P. The missing term in effective pair potentials. *J. Phys. Chem.* **91**, 6269–6271 (1987).

50. Darden, T., York, D. & Pedersen, L. Particle mesh ewald: an *n log(n)* method for ewald sums in large systems. *J. Chem. Phys.* **98**, 10089–10092 (1993).

51. Essmann, U. et al. A smooth particle mesh ewald method. *J. Chem. Phys.* **103**, 8577–8593 (1995).

52. Berendsen, H. J. C., Postma, J. P. M., van Gunsteren, W. F., DiNola, A. & Haak, J. R. Molecular dynamics with coupling to an external bath. *J. Chem. Phys.* **81**, 3684–3690 (1984).

53. Berger, O., Edholm, O. & Jähnig, F. Molecular dynamics simulations of a fluid bilayer of dipalmitoylphosphatidylcholine at full hydration, constant pressure, and constant temperature. *Biophys. J.* **72**, 2002–2013 (1997).

54. Oostenbrink, C., Villa, A., Mark, A. E. & Van Gunsteren, W. F. A biomolecular force field based on the free enthalpy of hydration and solvation: the GROMOS force-field parameter sets 53A5 and 53A6. *J. Comput. Chem.* **25**, 1656–1676 (2004).

55. Widom, B. Some topics in the theory of fluids. *J. Chem. Phys.* **39**, 2808–2812 (1963).

56. Frenkel, D. & Smit, B. *Understanding Molecular Simulation: From Algorithms to Applications* (Academic Press, 2001).

57. Schlaich, A., Kowalik, B., Kanduč, M., Schneck, E. & Netz, R. R. *Computational Trends in Solvation and Transport in Liquids, IAS Series* 155–185 (Forschungszentrum Juelich, 2015).

58. Waheed, Q. & Edholm, O. Undulation contributions to the area compressibility in lipid bilayer simulations. *Biophys. J.* **97**, 2754–2760 (2009).

## Acknowledgements

We thank Peter Fratzl for insightful comments and Institut Laue-Langevin for neutron beamtime. E.S. acknowledges financial support by the Max Planck Society and from an Emmy-Noether grant (SCHN 1396/1) of the German Research Foundation (DFG). R.R.N. acknowledges support of the DFG via SFB 1112.

## Author contributions

M.K., R.R.N. and E.S. designed research; M.K. and E.S. performed research; M.K., A.S. and A.H.V. developed analytical tools and carried out the analysis; J.J., E.M. and B.D. contributed experimental measurements; M.K., E.M., B.D., R.R.N. and E.S. wrote the paper.

## Additional information

*Supplementary Information* accompanies this paper at http://www.nature.com/naturecommunications.

**Competing interests:** The authors declare no competing financial interests.

**Reprints and permission** information is available online at http://npg.nature.com/reprintsandpermissions/

**How to cite this article:** Kanduč, M. et al. Tight cohesion between glycolipid membranes results from balanced water-headgroup interactions. *Nat. Commun.* **8**, 14899 (doi: 10.1038/ncomms14899 (2017)).

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_14.png" /></figure>

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

© The Author(s) 2017