We have investigated the interactions of Le^X and of Lac 2 glycolipids in simulations in which these glycolipids are embedded in POPC lipid membranes. Our Le^X and Lac 2 glycolipids have the same lipid tails as POPC and carbohydrate tips that are connected to these lipid tails by a glycerol linker group (see Figure 2). Standard carbohydrate force fields overestimate attractive carbohydrate–carbohydrate interactions, which leads to osmotic pressures for solutions of neutral carbohydrates that are systematically too low compared to experimental values (Lay et al., 2016; Sauter and Grafmüller, 2016). We have therefore used the G L Y C A M 06 OSMOr 14 TIP 5 P force field for the carbohydrate tips of our glycolipids, in combination with the standard AMBER Lipid14 force field (Dickson et al., 2014) for the glycolipid tails and the POPC lipids (Kav, 2019; Kav et al., 2020). In the G L Y C A M 06 OSMOr 14 TIP 5 P force field, the van der Waals parameters for carbohydrate–carbohydrate interactions of the standard force field GLYCAM06 have been reparametrized to correctly reproduce experimentally measured osmotic pressures (Sauter and Grafmüller, 2016). The G L Y C A M 06 OSMOr 14 TIP 5 P force field employs the TIP5P water model that leads to more reliable carbohydrate–carbohydrate interactions in GLYCAM06 than the standard TIP3P water model (Sauter and Grafmüller, 2015; Woods, 2018). Because simulations of AMBER Lipid14 POPC membranes in TIP5P water lead to an unreasonably small area per lipid, we have rescaled the Lennard–Jones interactions between the TIP5P water molecules and the POPC headgroup atoms to obtain the same area per lipid as in standard AMBER Lipid14 simulations with the TIP3P water model (Kav et al., 2020).

Our computational investigation is based on four different simulation systems. We have investigated the trans-interaction of a single pair of Le^X or Lac 2 glycolipids in a simulation system that consists of a lipid bilayer with one glycolipid anchored in each monolayer (Kav et al., 2020) (see Figure 1A). In this system, the two glycolipids in the different monolayers interact due to the periodic boundary conditions of the simulation box. We have varied the separation of the membrane monolayers by varying the number of water molecules in the simulation box. At each membrane separation, we have generated 10 trajectories with a length of 3 μs for the Le^X system and with a length of 1 μs for the Lac 2 system at a temperature of 30°C. The total simulation times at each membrane separation are thus 30 and 10 μs for the Le^X and Lac two systems, respectively. The membranes contain in each monolayer 35 POPC lipids besides the single glycolipid and have an area of 23.3 nm². Similarly, we have investigated the cis-interaction of a single pair of Le^X glycolipids in a simulation system in which two glycolipids are anchored in the same monolayer of the membrane (see Figure 1B). The monolayer with the two Le^X glycolipids contains 34 POPC lipids, and the other monolayer contains 36 POPC lipids. We have generated 25 trajectories with a length of 3.6 μs for this system at a temperature of 30°C. The total simulation time for this system is thus 90 µs.

In addition, we have investigated the interplay of trans- and cis-interactions in a system with 10 Le^X or 10 Lac 2 glycolipids in each monolayer of a membrane, besides 90 POPC lipids (see Figure 1C). The fraction of glycolipids in the membrane is thus 10 mol%, and the area of the membrane is 63.5 nm². In this system, the glycolipids cis-interact with glycolipids of the same monolayer and trans-interact with glycolipids in the other monolayer across the periodic boundaries of the simulation box. By varying the number of water molecules in the simulation box, we have varied the membrane separation (see Figure 3) and have generated 10 trajectories with a length of 1 μs at each separation. The values for the membrane separation in Figure 3 correspond to the separation from membrane midplane to membrane midplane and, thus, to the height of the simulation box. In the simulation systems of Figures 1A,C, the local membrane separation along the membranes is constant and equal to the simulation box height. The local separation of the membranes is not affected by the small membrane shape fluctuations in these systems because the two monolayers of the membrane are coupled in these shape fluctuations.

Finally, we have explored the role of membrane shape fluctuations in glycolipid-mediated adhesion in our largest simulation system, which consists of two membranes with 90 Le^X glycolipids and 810 POPC lipids in each of the four monolayers of the membranes (see Figure 1D). In this system, each membrane has an area of 582 nm², and the two aqueous compartments between the membranes both contain 70,776 water molecules so that the average membrane separation is close to the average separation measured in neutron scattering experiments with stacks of membranes that contain 10 mol% of Le^X glycolipids (Schneck et al., 2011). The local separation between the two membrane interfaces varies in this system because of the shape fluctuations of the two membranes (see Figure 1D). We have generated 10 trajectories of this system with lengths that vary between 0.95 and 1.12 μs and sum up to a cumulative simulation time of 10.43 µs

We have generated initial structures of the simulation systems in Figures 1A–C by first building POPC lipid membranes using the CHARMM-GUI program (Jo et al., 2008) and subsequent replacement of lipids by Le^X or Lac 2 glycolipids, which have the same lipids tails as POPC. After initial minimization and equilibration following standard procedures (Dickson et al., 2014; Kav et al., 2020), we have produced all simulation trajectories in AMBER 16 GPU (Salomon-Ferrer et al., 2013b; Le Grand et al., 2013) at a simulation temperature of 30°C using a Langevin thermostat (Salomon-Ferrer et al., 2013a) with a collision frequency of 5.0 ps⁻¹. We have employed a semi-isotropic pressure coupling with a pressure of 1 bar in all directions, which corresponds to a membrane tension of zero, and the Berendsen barostat (Berendsen et al., 1984) with a relaxation time of τ = 3 ps for the pressure regulation because of the stability of the semi-isotropic pressure coupling in AMBER 16 GPU with this barostat. We have constrained the bond lengths for hydrogen atoms using the SHAKE algorithm (Ryckaert et al., 1977; Miyamoto and Kollman, 1992) and have set the MD integration timestep to 2 fs. The nonbonded interactions were calculated using the Particle Mesh Ewald (PME) algorithm (Darden et al., 1993; Essmann et al., 1995) for a cutoff length of 1.0 nm.

To generate an initial structure for the large simulation system of Figure 1D, we have used the initial structure for the simulation system of Figure 1C, which contains 10 Le^X glycolipids and 90 lipids in each monolayer. We have first replicated this initial structure 9 times to generate the first membrane with 90 Le^X glycolipids and 810 lipids in each monolayer and have then duplicated and translated the membrane to obtain the second membrane. In order to increase the computational efficiency, we have applied hydrogen mass repartitioning (Hopkins et al., 2015), which allowed us to increase the MD integration time step to 3 fs, in addition to the equilibration and simulation procedures described above for the smaller systems.

We have identified interaction events between the carbohydrate tips of the glycolipids along the simulation trajectories as consecutive stretches of simulation frames at intervals of 0.1 ns with nonzero contacts of the tips (Kav et al., 2020). Here, contacts are defined as contacts between non-hydrogen atoms of the two carbohydrate tips within a distance of less than 0.45 nm. The interaction events of the carbohydrate tips can be characterized by their lifetime and by the maximum number of contacts of the events. In an attempt to distinguish between collisions and binding events, only interaction events with a maximum number of contacts that is larger or equal to a cutoff number n _c are considered as binding events. This distinction based on a cutoff number n _c of contacts is somewhat arbitrary because of the fuzzy interactions of the carbohydrates (Kav et al., 2020), which exhibit a large variety of diverse, bound conformations in our simulations, rather than a single binding conformation (see Figure 4). The binding constants determined from our simulations therefore depend—to some extent—on the cutoff number n _c . All binding constant values reported here have been calculated for n _c = 5. Binding constant values for n _c = 10 are typically about 10–15% smaller than the values obtained for n _c = 5 (Kav et al., 2020).

We have determined the trans-binding constant K _trans of the two Le^X or two Lac 2 glycolipids embedded in the different membrane monolayers of the simulation system shown in Figure 1A as K _trans = AP _b /(1 − P _b ), where P _b is the probability that the two glycolipid tips are bound, and A is the membrane area. The binding probability P _b is simply the fraction of simulation frames belonging to binding events. We have calculated the trans-binding probability of the Le^X glycolipids from 1,490, 1,609, 588, and 141 binding events observed on our trajectories at the membrane separations 5.5, 6.0, 6.5, and 7.0 nm, respectively. For two Lac 2 glycolipids, we have obtained 609, 413, 183, and 34 such binding events on the trajectories at the corresponding membrane separations. To ensure independence from the initial conformation of the trajectories, we have discarded the first 10% of each trajectory in our calculations of the binding probability P _b of the two molecules. Similarly, we have determined the cis-binding constant K _cis of the two Le^X glycolipids embedded in the same membrane monolayer of the simulation system of Figure 1B as K _cis = AP _b /(1 − P _b ), where P _b is the binding probability obtained from 4,308 binding events on all trajectories, after discarding the first 20% of each trajectory to ensure independence from the initial conformation.

In the simulation system of Figure 1C, a Le^X or Lac 2 glycolipid can be involved in several trans- and/or cis-interactions in a simulation frame. In this system, we have determined the effective trans-binding constant as K trans = A n trans / ( 10 − n trans ) 2 where n _trans is the average number of glycolipids in a monolayer engaged in trans-interactions. Because each of the monolayers contains 10 glycolipids, the average number of glycolipids that are not engaged in trans-interaction is 10 − n _trans. For all our simulations, errors have been calculated as the error of the mean of values obtained for the independent trajectories.

Experiments were carried out with DPPC membranes doped with 10 mol% glycolipids, whose synthesis has been described elsewhere (Schneider et al., 2001; Schneider et al., 2003; Tanaka et al., 2003; Tanaka et al., 2004). In the glycolipids used in the experiments, the alkyl chains were connected to the glycerol via ether bonds, in contrast to the glycolipids in the simulations, where the alkyl chains are connected via ester bonds. The Lac 1 glycolipid (Schneider et al., 2001, Schneider et al., 2003; Tanaka et al., 2003; Schneck et al., 2008), Lac 2 glycolipid (Schneider et al., 2001, Schneider et al., 2003; Tanaka et al., 2003), and Gentiobiose glycolipid (Tanaka et al., 2004; Schneck et al., 2008) consist of two saturated hexadecyl chains and of carbohydrate tips connected to the chains via a glycerol junction (see Supplementary Figure S1). The Lac 1 and Gentiobiose glycolipids have the disaccharides lactose and gentiobiose as carbohydrate tips, while the Lac 2 glycolipid has a tetrasaccharide tip composed of two lactose units. Neutron diffraction contrast was optimized by using chain-deuterated DPPC and light water (H₂O) in combination with the Lac 1 and Gentiobiose glycolipids, which were available in their chain-deuterated forms. Consequently, chain-hydrogenous DPPC and heavy water (D₂O) were used in combination with the Lac 2 glycolipid, which was available in its chain-hydrogenous form. All lipids were initially dissolved in 7:3 mixtures (v/v) of chloroform and methanol at a concentration of 1 mg/ml. Subsequently, mixed solutions with the desired molar ratio were prepared. Of these solutions, 1–2 ml were deposited onto planar Si[100]-substrates with native oxide (Si-Mat, Landsberg/Lech, Germany), which had previously been cut into a rectangular shape (65 mm × 25 mm) and cleaned using a modified RCA method (Kern and Puotinen, 1970). Due to their amphiphilic nature, the lipid mixtures form aligned membrane multilayers on the planar surfaces. To remove residual solvent, the coated wafers were stored at 70°C for 3 h and, subsequently, in a vacuum chamber overnight. At least two heating/cooling cycles between 20 and 80°C were performed at a high relative humidity > 95 %.

Neutron diffraction experiments were performed on the high-resolution diffractometer D16 at the Institut Laue-Langevin (ILL, Grenoble, France). The incident beam with a wavelength of λ = 0.474 nm (Δλ/λ ≈ 1%) reached the sample plane through the aluminum windows of the sample chamber, with an adjustable angle of incidence Ω. Scattering occurs into various directions at angles Γ with respect to the incident beam. For each Ω, the Γ-dependent intensity is recorded using a position-sensitive ³He detector with 128 × 128 channels and a spatial resolution of 2 mm. By rotating the sample stage and, thus, by stepwise variation of Ω, two-dimensional maps of the intensity as a function of Γ and Ω were recorded (see images in Supplementary Figure S2). During this procedure, the intensity was normalized to the intensity of the incident beam (via an in-beam monitor), the pixel sensitivity and solid angle, and the illuminated sample area. Bragg peaks associated with the lamellar period D of the membrane multilayers are found where the specular condition (Γ = 2Ω) coincides with the Bragg condition (Γ = 2 arcsin(nλ/(2D)), see plots in Supplementary Figure S2), where n ∈ N . The average membrane separation is identical to the lamellar period and obtained by solving the Bragg condition for D.

For the diffraction experiments, a liquid cell designed for solid-supported membrane multilayers (Schneck et al., 2008, Schneck et al., 2009) was used: two planar Si substrates, one of them coated with the membrane multilayers, were assembled into a sandwich-like configuration with small glass pieces (thickness: 0.10 mm) as spacers between them. The space between the two wafers was then filled with H₂O- or D₂O-based aqueous solutions containing 100 mM NaCl and 5 mM Hepes (Fluka, Taufkirchen, Germany) and, optionally, 5 mM CaCl₂. During the diffraction experiments, the measurement cell was placed in a climate chamber at controlled temperature and high relative humidity ( > 95 %) to minimize water evaporation. Temperature was kept at 60°C such that the bilayers were in the fluid L _α -phase in all cases.

Figure 5A illustrates the trans-binding constant K _trans of Le^X glycolipids and of Lac 2 glycolipids obtained in the simulation system of Figure 1A with a single trans-pair of glycolipids and in the system of Figure 1C with 10 mol% glycolipids. In both systems, the glycolipids in the different monolayers trans-interact due to the periodic boundary conditions of the simulation box. In the system of Figure 1C, glycolipids in the same monolayer cis-interact in addition to the trans-interactions with glycolipids in the apposing monolayer. We have determined K _trans at different membrane separations by varying the number of water molecules in the simulation box (see Figure 5B). The values for the membrane separation l given in Figure 5 correspond to the separation from membrane midplane to membrane midplane and, thus, to the height of the simulation box. At the smallest separation (l = 5.5–5.6 nm) considered in the different simulation systems, the number of water molecules per lipid in the simulation box ranges from 16 to 19, depending on the system. Simulation conformations of the membranes with 10 mol% Le^X glycolipids or Lac 2 glycolipids at the different separations are shown in Figure 3.

The trans-binding constants K _trans of the glycolipids strongly decrease with increasing membrane separation. This decrease in trans-interactions can be understood from a decrease in the overlap of the carbohydrate tips of the glyoclipids. Figure 4 illustrates 50 randomly selected trans-complexes of the Le^X glycolipid tips with at least 10 contacts of non-hydrogen atoms at each of the membrane separations l = 5.5, 6.0, 6.5, and 7.0 nm. The carbohydrate tip of the lower Le^X glycolipid is aligned in the 50 complexes and represented in blue colors, while the carbohydrate tip of the upper glycolipid is represented in red/yellow colors. The clouds of red/yellow carbohydrates illustrate that the interactions of the glycolipids are fuzzy, that is, the carbohydrate tips exhibit a large variety of bound conformations in our simulations. The overlap of the cloud of the upper, red/yellow carbohydrates with the lower, blue carbohydrate decreases with increasing membrane separation. At the membrane separation l = 5.5 nm, the Le^X glycolipids interact via their entire carbohydrate tips. At l = 6.0 nm, the interactions are limited to the branched trisaccharide of the glycolipid tip, and at l = 6.5 and 7.0 nm, the interactions are further restricted to the galactose and fucose monosaccharides at the branched end of the Le^X glycolipid.

A comparison of the K _trans values from the two simulation systems indicates that cis-interactions tend to strengthen the trans-interactions. At a given membrane separation, the K _trans values obtained in the system with 10 mol% glycolipids, in which the glycolipids can also cis-interact, are slightly larger than the K _trans values obtained in the system with a single trans-pair of glycolipids. For Le^X glycolipids, we have determined the cis-binding constant K _cis = 7.9 ± 1.2 nm² in the simulation system of Figure 1C with a single cis-pair. This value of K _cis is comparable to and even slightly larger than the K _trans values at the membrane separation 5.5 nm, at which the Le^X glycolipids can trans-interact with their entire carbohydrate tips according to Figure 4. The comparable magnitude of cis- and trans-interactions at this membrane separation is also reflected by the average number of Le^X glyolipids engaged in trans- and cis-interactions in the membrane system with 10 mol% of glycolipids. We find that, on average, about four out of the 10 Le^X glyolipids of a monolayer are engaged in trans-interactions, while about three of the 10 Le^X glyolipids are engaged in cis-interactions at the smallest membrane separation considered in our simulations. Because of the fuzzy binding, trans- and cis-interactions of the glycolipids are not mutually exclusive, and trimer and higher multimers in different trans- and cis-combinations frequently occur along the simulation trajectories.

In addition to Le^X and Lac 2 glycolipids, we have also determined the trans-binding constant K _trans of Lac 1 glycolipids with a small disaccharide galactose–glucose tip in the system of Figure 1C at the membrane separations l = 5.3, 5.6, 6.0, and 6.5 nm. The carbohydrate tip of the Lac 1 glycolipid thus has half the size of the tetrasaccharide tip of the Lac 2 glycolipid. For these Lac 1 glycolipids, we only obtain noticeable trans-interactions at the smallest membrane separation 5.3 nm with K _trans = 0.90 ± 0.07 nm². This value of K _trans is comparable to the K _trans value of Lac 2 glycolipids at the separation 6.5 nm in the same system with 10 mol% of glycolipids (see Figure 5A). At the membrane separation 5.6 nm, the K _trans value of the Lac 1 glycolipids is already strongly reduced to 0.06 ± 0.01 nm² and vanishingly small at the separations 6.0 and 6.5 nm.

Figure 6 illustrates neutron scattering results for the average separation in stacks of DPPC membranes that contain 10 mol% of Lac 1, Gentiobiose, Lac 2, or Le^X glycolipids. The neutron scattering experiments have been performed at a temperature of 60°C at which the DPPC membrane is fluid. We find that the average separation between adjacent membranes in the stack strongly depends on the glycolipid embedded in the membranes. In the absence of Ca²⁺ ions that associate with PC lipid headgroups and induce an electrostatic repulsion of the membranes (Lis et al., 1981; Altenbach and Seelig, 1984), the average membrane separation of DPPC membranes that contain 10 mol% of Lac 1 or Gentiobiose lipid with small disaccharide tips is close to the average separation of pure DPPC membranes. The average separations of membranes with 10 mol% of Le^X and Lac 2 glycolipids, in contrast, are considerably larger than the average separation of pure DPPC membranes. In the presence of 5 mM Ca²⁺ ions, the average separation of pure DPPC membranes strongly increases by 2.4 nm, whereas the average separation of membranes with 10 mol% of Lac 1, Gentiobiose, and Le^X glycolipids only increases by about 0.3–0.4 nm, compared to the separation in the absence of Ca²⁺. Interestingly, the average separation of membranes with 10 mol% of Lac 2 glycolipids increases to a value close to the separation of pure DPPC membranes in the presence of 5 mM Ca²⁺.

For Le^X glycolipids, the average membrane separation from neutron scattering has been previously reported at concentrations of 0 mol% (pure DPPC), 2 mol%, 5 mol%, 10 mol%, and 25 mol% both in the absence of Ca²⁺ and the presence of 1, 2, and 5 mM Ca²⁺ (Schneck et al., 2011). At the different Ca²⁺ concentrations, the average separation of the membranes attains values that are rather similar at large concentrations of the Le^X glycolipids, which indicates that these average separations result from stable membrane adhesion mediated by the Le^X glycolipids. The average membrane separations are significantly larger than the separations at which Le^X glycolipids trans-interact in our simulations. Membrane shape fluctuations therefore seem to be important to bring the glycolipids in the apposing membrane surfaces into local contact.

We have investigated the role of membrane shape fluctuations in our largest simulation system that consists of two membranes with 10 mol% of Le^X glycolipids and an area of 24.1 × 24.1 nm². In this system, we have adjusted the average separation from membrane midplane to midplane to 7.7 nm, the average membrane separation measured in the neutron scattering experiments. Using the same separation value is justified because fluid DPPC membranes (Lis et al., 1982) and the POPC membranes simulated here have virtually the same thickness when the water layer thickness is defined consistently through the water volume and the area per lipid, as described in the study by Kanduc et al. (2017). In our simulations, the two aqueous compartments between the membranes contain the same number of water molecules, and the simulation box height of 15.4 nm is twice as large as the average membrane separation because we have two membranes. The membrane shape fluctuations can be quantified by the relative roughness of the two membranes, which is the standard deviation ξ ⊥ = ⟨ ( l − l ̄ ) 2 ⟩ of the local separation l of the membranes from the average separation l ̄ . To calculate the roughness ξ _⊥ of the membranes in a simulation conformation, we have divided the x-y-plane of our simulation box, which is, on average, parallel to the membranes, into 16 × 16 quadratic patches. This discretization of the x-y-plane leads to 16 × 16 membrane patches in each of the two membranes. The membrane patches have an area of about 1.5 × 1.5 nm² and contain, on average, 3.5 POPC lipids or glycolipids in each monolayer. For each patch in each of the two membranes, we have determined the z-coordinate of the center of mass of the lipid tails in the same monolayer and have calculated the z-position of the membrane midplane as the average value of the z-coordinates for the monolayers of the membrane patch. From the two z-positions of the two apposing membrane patches with the same x-y-position, we have obtained two values for the local separation l between these membrane patches, which add up to the simulation box height.

Figure 7A illustrates the membrane roughness along four of the 10 simulation trajectories for this system. Depending on the trajectory, the membrane roughness varies between about 0.3 and 1 nm, after an initial relaxation of about 0.15 μs, during which the roughness increases on all trajectories. The last conformations of these four trajectories are shown in Figure 8. The roughness of the conformation at the top right of Figure 8, which corresponds to the yellow trajectory in Figure 7A, is visibly larger than the roughness of the other three conformations and leads to trans-interactions of Le^X glycolipids in the lower membrane interface. The interplay of membrane roughness and trans-interactions is also supported by the correlation between the mean membrane roughness and the mean number of trans-bonds of the 10 trajectories after the initial relaxation of about 0.15 µs (see Figure 7B).

As expected from Figure 5, trans-bonds of Le^X glycolipids along the trajectories only occur at local membrane separations l that are smaller than the average separation l ̄ = 7.7 nm of the membranes (see Figure 9A). The distribution of the local membrane separations for membrane patches with trans-bound Le^X glycolipids is centered around 6.7 nm, while the distribution for all membrane patches is centered at the average membrane separation of 7.7 nm of the membranes. The standard deviation of the distribution P(l) is the membrane roughness ξ _⊥ = 0.53 nm averaged over all trajectories. This average membrane roughness along our trajectories is somewhat smaller than the membrane roughness ξ _⊥ = 0.73 ± 0.03 nm obtained from neutron scattering experiments (Schneck et al., 2011; Kav et al., 2020), likely because the membrane area in our large simulation system is still somewhat too small to allow all relevant fluctuation modes of the membranes. The lateral correlation function of the local separation in Figure 9B indicates remaining correlations between membrane patches with the maximum distance of 12 nm along the x- and y-axis of the simulation box, which lead to negative values of the correlation function at this maximum distance.