Free energies of membrane stalk formation from a lipidomics perspective

Many biological membranes are asymmetric and exhibit complex lipid composition, comprising hundreds of distinct chemical species. Identifying the biological function and advantage of this complexity is a central goal of membrane biology. Here, we study how membrane complexity controls the energetics of the first steps of membrane fusions, that is, the formation of a stalk. We first present a computationally efficient method for simulating thermodynamically reversible pathways of stalk formation at coarse-grained resolution. The method reveals that the inner leaflet of a typical plasma membrane is far more fusogenic than the outer leaflet, which is likely an adaptation to evolutionary pressure. To rationalize these findings by the distinct lipid compositions, we computed ~200 free energies of stalk formation in membranes with different lipid head groups, tail lengths, tail unsaturations, and sterol content. In summary, the simulations reveal a drastic influence of the lipid composition on stalk formation and a comprehensive fusogenicity map of many biologically relevant lipid classes.

, an older MARTINI POPC model with a 5-bead oleoyl tail was used for both simulations, longer than the four-bead oleoyl model used for all other simulations of this study. The systems contained 128 POPC lipids per bilayer and either 200 (PChd200) or 220 (PChd220) water beads in the proximal compartment, corresponding to 1.56 or 1.72 water beads per lipid. Evidently, the PMFs along ξ ch suggest lower free energies for the stalk as compared to the MFEPs along the density-based order parameter used in Ref. 1. The difference of ∼30 kJ/mol may be rationalized by the different definitions of the end states. In our method, the stalk end state (ξ ch ≈ 1) is defined by the presence of a hydrophobic connection between the two membranes, but all possible shapes, radii, and lateral positions of the connection are included in the state with ξ ch ≈ 1. In Ref. 1, the stalk state is defined with a specific 3D density of the hydrophobic beads, which may allow fewer conformational states than the stalk definition adopted by us. 2 Source data are provided as a Source Data file.   Mass densities at lipid stalk (set 1) for systems POPC:Cholesterol 80:20, POPC:LPC 80:20, POPC:PCN 80:20, POPC:PCA 80:20. Densities were computed from the umbrella window restrained to ξ ch = 1, omitting the first 50 ns for equilibration. The densities ρ(r, z) were computed as function of the lateral distance r and normal distance z from the center of the stalk, defined as the center of the cylinder used to define ξ ch (see Methods). ρ(r, z) were copied to negative r-values purely for visualization purposes.

Simulation setup and parameters
The double-membrane systems were set up with a multi-step protocol that ensures the requested degree of hydration between the two membranes, irrespective of the lipid composition. First, a single membrane system was set up with the Insane software. 3 The membrane was hydrated with the requested number of water beads and neutralized with sodium beads, as needed. The energy of the system was minimized, and the membrane was equilibrated for 20 ns to fully relax the box dimensions.
Next, two copies of the membrane were stacked on top of each other, (i) leading to the requested degree of hydration between the proximal leaflets of the double-membrane and (ii) to balanced electric charges between the two water compartments, thereby avoiding a transmembrane electric potential. To fully hydrate the distal leaflets, the box was enlarged along the z direction and additional water was added to the distal water compartment, until the distal leaflets were hydrated with 10 water beads per lipid. The double membrane was equilibrated for another 20 ns. The simulation setup was automated with Bash scripts.
In long simulations of the double-membrane system, we Visual inspection of the simulation, also after long simulation times, did not reveal any indication of phase separation or lipid demixing (see Fig. 5).

Number of lipids per leaflet and membrane size
To exclude finite-size effects from the periodic boundaries, we computed PMFs of stalk formation for three lipid compositions (pure DOPC, pure DIPC, or a 70:30 mixtures of DOPC with cholesterol) at various membrane sizes (Fig. 1). For the pure-lipid systems, only marginal finite-size effects were evident from the PMFs for very small systems (Fig.   1D, black and orange). Hence, systems without cholesterol were built with 64 lipids per monolayer throughout this study. In contrast, for cholesterol-containing membranes, the stalk free energies were increased by up to 30 kJ/mol for very small box sizes relative to larger boxes, likely because cholesterol renders the membranes more rigid (Fig. 1D, blue).
Hence, cholesterol-containing membranes were built with simulation box sizes of approx.
7.5 nm in the x-and y-dimensions (visualized in Fig. 5).

Reaction coordinate for stalk formation
PMFs were computed along the "chain coordinate" ξ ch , which was originally introduced to drive pore formation in membranes. 8,9 ξ ch quantifies the connectivity of two compartments of specific atoms. In this study, ξ ch quantified the degree of connectivity between the hydrophobic cores of two membranes. The shape, radius, lateral position, or lipid composition of the hydrophobic connection were not controlled ξ ch but decided by the force field.
ξ ch was defined by a cylinder of radius 1.2 nm that spans the two hydrophobic regions of the two membranes and the proximal water compartment (Main Article Figure 1A). The cylinder is decomposed into slices with a thickness of 1Å. Then, ξ ch is approximately given by the fraction of slices that are filled by lipid tail beads: Here, N s is the number of cylinder slices, n s ≥ 0). Hence, upon pulling the system along ξ ch , the slices are filled one-by-one, thereby gradually forming a hydrophobic connection between the two hydrophobic membrane cores, as required for stalk formation. As bead contributing to ξ ch , hydrophobic lipid tail beads as well as hydrophobic beads of cholesterol were used.
Critically, the radius of the cylinder of 1.2 nm does not control the radius of the stalk.
Instead, the cylinder is merely used to ensure the locality of the hydrophobic protrusions in the membrane plane. If the cylinder radius would be too large, two laterally displaced hydrophobic protrusions, one from the upper and one from the lower membrane, could be misinterpreted as a continuous hydrophobic connection. Such problems would further allow the membrane to evade the energetically unfavorable transition state of stalk formation, which could lead to undesired hysteresis effects. 10 Figure 17 shows that the PMFs depend only marginally on the choice of the cylinder radius in the range between 1.2 nm and 1.6 nm. The lateral position of the cylinder was not fixed but dynamically defined to allow the cylinder to "follow" the stalk as the stalk travels parallel to the membrane plane. This property excludes that the system moves along ξ ch by shifting the stalk laterally out of the cylinder, which would again lead to undesired hysteresis effects. 11 For implementation details, we refer to previous work. 8 To visualize the stalk in molecular graphics, the stalk was translated to the box center purely for illustration purposes (Figs. 1, 5, 8 9, 10).
To render ξ ch differentiable, the function δ ζ was defined with a differentiable switch function that approximates and indicator function: 8 Here, the parameter ζ indicates the fraction to which the slice is filled upon the addition of the first apolar bead into the slice. We used ζ = 0.75 in this study. The parameters b and c are taken as b = ζ/(1 − ζ) and c = (1 − ζ)e b , leading to a continuous and differentiable switch function. Likewise, the number of beads n (t) s in slice s was defined with a differentiable indicator function: Here, the sum is taken over all N b apolar beads and r j denotes the Cartesian coordinates of bead j. f is a three-dimensional indicator function that takes unity inside the volume of slice s, and f smoothly switches to zero at the slice boundaries.

Definition of the differentiable indicator function
The indicator function f (r j ) was formulated as a product of an axial (along z direction) and radial indicator function (in lateral direction), as follows: 8 Here, r j = (x j , y j , z j ) are the Cartesian coordinates of atom j. X cyl and Y cyl is the position of the cylinder axis in the membrane x-y-plane, and R cyl = 1.2 nm is the radius of the cylinder. The radial and axial indicator functions are defined with the help the following differen-tiable step function θ(x; h), where 2h is the width of the switch region: Hence, θ(x; h) is zero at |x| > 1 + h and unity at |x| < 1 − h with continuous switches in the regions 1 − h ≤ |x| ≤ 1 + h. In this work, we used h = 1/4. The axial and the radial indicator functions were defined as where r j = [(x j − X cyl ) 2 + (y j − Y cyl ) 2 ] 1/2 denotes the lateral distance of atom j from the cylinder axis.
The z s coordinates of the cylinder slices are defined relative to the center of mass Z prox of the proximal water beads, keeping the cylinder centered between the two membranes. The position of the cylinder axis (X cyl , Y cyl ) is defined via a weighted center-of-mass calculation of hydrophobic tail beads in the lateral layer spanned by the cylinder height, i.e., by tail beads whose z coordinate fulfill Z prox − N s d/2 ≤ z ≤ Z prox + N s d/2. This definition ensures that that cylinder follows the stalk as it travels in the lateral plane. More details as well as the inner derivatives required for computing the forces derived from restraints along ξ ch were described previously. 8

Umbrella sampling simulations of stalk formation
PMFs were computed using umbrella sampling (US). Initial frames for US were taken from constant-velocity pulling simulations, in which the systems were pulled from ξ ch = 0.1 to ξ ch = 1 within 200 ns, using a force constant of 3000 kJ/mol. Visual inspection of the simulations showed that pulling ξ ch led to gradual stalk formation in all membrane system.
19 umbrella windows were used with reference positions between 0.1 and 1 in steps of 0.05.
The force constant was set to 3000 kJ/mol. Each window was simulated for 200 ns, where the first 50 ns were omitted for equilibration. An integration time step of 20 fs was used.
All other parameters were chosen as described above. The PMFs were computed with the weighted histogram analysis method (WHAM), 12 as implemented in the gmx wham module of Gromacs. 13 The free energy of stalk formation ∆G stalk was defined as the value of the PMF at ξ ch = 0.95.

Statistical errors of PMFs and ∆G stalk
Statistical errors were estimated with the Bayesian bootstrap of complete histograms. 13 Accordingly, in each round of bootstrapping, random weights were assigned to all histograms, and the randomly weighted histograms were used to compute a bootstrapped 'synthetic' PMF. Each bootstrapped PMF was defined to zero at ξ ch = 0.

Absence of hysteresis between stalk opening and closing pathways
As a test for the validity of the reaction coordinate, we computed the PMFs along stalkopening and stalk-closing pathways (Fig. 18). We obtained nearly identical PMFs along opening and closing pathways, confirming the absence of undesired hysteresis.

Unbiased simulations of stalk formation and closure
To test whether the PMF along ξ ch reflects the true free energy difference between the flat membrane and the stalk, and to obtain the rates of stalk formation and closure, we carried out unbiased simulations. Here, we used the beta-3.2 release of MARTINI 3.0. We simulated a double-membrane of pure POPC with 230 and 1920 water beads in the proximal and distal water compartments, respectively, for which the PMF suggested a free energy difference between stalk and flat membrane of ∆G stalk ≈ 0 (Main Article Fig. 2). Four replicas of 200 µs each were simulated, which carried out 8 transitions of stalk formation and 7 transitions of stalk closure corresponding to rates of k stalk = 16 ms −1 and k closure = 23 ms −1 , respectively. Hence, the free simulations suggest a free energy of stalk formation of ∆G stalk = −k B T ln(k stalk /k closure ) = 0.9 kJ/mol, in excellent agreement with the PMF.

Density calculations
The mass density around the stalk was computed with an in-house modification of the

PMF calculations of dehydration of the proximal leaflet
PMFs for dehydration (Fig. 12) were computed similar to the work by Smirnova et al. 1,14 The double-membrane system was setup and equilibrated as described above using 26 water beads per nm 2 and 64 lipid per leaflet. The membranes were composed of pure DOPC or of DOPC/cholesterol mixtures with molar ratios of 80:20 or 60:40. The reaction coordinate ξ dehyd was taken as the center-of-mass distance in z direction between the phosphate beads of the two proximal leaflets (Fig. 12, black arrows). To allow the equilibration of water beads between the two water reservoirs, we introduced a small hole in both leaflets. The hole allowed occasional permeation of water beads, thereby keeping the water beads of the proximal compartment at approximately constant chemical potential. 1, 14 We implemented the hole using a repulsive cylindrical flat-bottomed potential acting on all lipid beads with a radius of 0.5 nm and a force constant of 1000 kJ mol −1 nm −2 .
The proximal compartment was dehydrated using constant-velocity pulling along ξ dehyd from ∼3.0 nm to 0.5 nm over 1 µs of simulation with a force constant of 4000 kJ mol −1 nm −2 .
The PMFs were computed using umbrella sampling and WHAM 12 with 56 equispaced um-brella windows at 0.05 nm distance and force constant of 2000 kJ mol −1 nm −2 . Starting frames were taken from the constant-velocity pulling simulation. Each window was simulated for 3 µs, where the first 1.5 µs were omitted from analysis for equilibration. All simulation parameters were chosen descried above. Statistical errors were computed by bootstrapping complete histograms as described above. 13

Supplementary Discussion
Attempt frequency of stalk formation coincides with the frequency of headgroup rearrangements along the membrane normal According to transition state theory, the rate of barrier crossing is given by where ν is the attempt frequency and ∆G ‡ = 25 kJ/mol the barrier height in the PMF  Fig. 19).
Taken together, the conformational sampling of POPC lipids studied here leads to both lateral and normal displacements on similar time scales of few nanoseconds. The lateral displacements manifest in lateral diffusion, whereas the normal displacements may be inter-preted as attempts for stalk formation with a success rate of e −∆G ‡ /k B T . It is important to note that, because the Martini model leads to a smoothed energy landscape, all dynamics discussed here (for lateral and normal displacements as well as for stalk formation) are likely accelerated relative to atomistic simulations or to experimental conditions. 15 Supplementary References