Free energies of stalk formation in the lipidomics era

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 near-atomic resolution. The new 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.


Introduction
Eukaryotic cellular membranes contain more than ten lipid classes, while each class comprises hundreds of different chemical species. 1,2 The complexity of membranes is further increased by the membrane asymmetry, that is, by distinct lipid compositions in the two leaflets. In the plasma membrane of mammals, sphingolipids are typically enriched in the outer leaflet of the plasma membrane, whereas phosphatidylethanolamine (PE) and phosphatidylserine (PS) are enriched in the inner leaflet. 3 Apart from different head groups, lipid species differ by the length and unsaturation of the fatty acid tails. Recent lipidomics studies showed that polyunsaturated lipids are strongly enriched in the inner leaflet. 4 Understanding why biological cells synthesize and maintain this complex lipid repertoire, that is, defining the biological function and advantage of specific lipid compositions, remains a central goal of membrane biophysics.
Here, we study how complex lipid compositions control the early stages of membrane fusion. Fusion is critical for many processes that involve the transport of cargos across membranes such as exocytosis, neurotransmission, infection of cells by enveloped viruses, fertilization, and intracellular transport. [5][6][7][8] Membrane fusion occurs with the help of fusion proteins because fusion requires the two opposing lipid membranes to overcome a large dehydration barrier. The process starts with two membranes in the lamellar phase at the equilibrium distance of 2-3 nm, followed by bridging of two opposing membranes by lipid acyl chains to establish a point-like protrusion. The lipids of the two contacting, proximal leaflets mix to form the initial transient hemifusion stalk. The stalk expands to allow the formation and expansion of a fusion pore. [9][10][11] Because intermediate structures along the fusion pathway involve highly curved membranes, the intrinsic curvature of lipids influences the kinetics of fusion. For instance, when added to the proximal leaflets, inverted cone-shaped lysophosphatidylcholine (LPC) inhibits stalk formation because it is incompatible with the large negative curvature at the stalk rim. Inversely, cone-shaped unsaturated PE lipids or diacylglycerol promote stalk forma-tion. 10,12,13 Here, the unsaturated fatty acids render the lipids more cone-shaped because the double bonds increase the tail disorder and, hence, the effective volume of the lipid tails. 13 However, in addition to such effect by the geometric lipid shape, the increased abundance of double bonds increases the conformational flexibility of the lipids, thereby allowing the membrane to adapt more easily to curvature, 14 which may further promote fusion.
Theoretical and computational approaches have established possible fusion pathways and the underlying free energy landscapes, however, only few studies focused on the role of the lipid composition during fusion. Early studies employed continuum descriptions, which model the role of lipids in terms of the effective spontaneous curvature and membrane rigidity, 15 as well as minimal coarse-grained lipids models in conjunction with Monte-Carlo, fieldtheoretic, or Brownian dynamics methods. [16][17][18] Complementary, to account for the chemical specificity of lipid-lipid interactions, fusion was studied with atomistic and coarse-grained (CG) molecular dynamics (MD) simulations, 19-21 discussed in several excellent reviews. [22][23][24] PE lipids were found to enhance fusion during simulations as expected from their negative intrinsic curvature, [25][26][27] however the effects of other head groups, sterols, or degrees of unsaturation have hardly been considered despite their abundance in biological membranes.
Computationally efficient free energy calculations of fusion require the definition of one or several reaction coordinates (or order parameters) along the fusion pathway; however, finding good coordinates for complex transitions is non-trivial. With poor coordinates, hysteresis problems emerge, barriers may be integrated out, and simulations proceed along non-reversible pathways. To avoid such problems, a recent elegant MD study used the string method to optimize the minimum free energy pathway of stalk formation, parametrized by the three-dimensional density of the apolar lipid tail beads. 27 However, because the string method is computationally demanding, it becomes prohibitive for high-throughput studies. In this study, we introduce a reaction coordinate for stalk formation that allows computationally highly efficient and thermodynamically reversible free energy calculations of stalk formation using common umbrella sampling simulations. The potentials of mean force (PMFs) along the coordinate, computed with Martini coarse-grained models, show that the inner leaflet of a common plasma membrane is far more fusogenic than the outer leaflet. To rationalize these finding, we screen the fusogenicity of lipids by varying the hydration level between the two membranes, the lipid headgroup size and charge, acyl chain length, and unsaturation levels. In addition to phosphatidylcholine (PC), PE, PS, and phosphatidylglycerol (PG) lipids, we also considered cholesterol, lysolipids, fatty acids in protonated and deprotonated forms, phosphatic acid, ceramide, diacylglycerol, and sphingomyelin. The calculations provide a comprehensive view on the role of lipid complexity during the first steps of membrane fusion.

Highly efficient free energy calculations of stalk formation
Potential of mean force (PMF) calculations may provide the free energy along reaction path in a computationally efficient manner. However, they require the choice of a good reaction coordinate (or order parameter); otherwise, undesired hysteresis effects may emerge, or energy barriers may be integrated out. We propose using a reaction coordinate for stalk formation that was originally introduced to study the formation of aqueous pores over membranes. 28,29 The coordinate, named "chain coordinate" ξ ch , quantifies the degree of connectivity between two compartments. In brief, ξ ch is defined with the help of cylinder that spans the head group and water regions of the proximal compartment (Fig. S1). The cylinder is decomposed into N s slices, and ξ ch is defined as the fraction of slices that are filled by apolar lipid tail beads. By pulling the system along the coordinate, the slices are filled one-by-one, thereby forming a continuous apolar connection between the two fusing membranes, as required for stalk formation ( Fig. 1A-C). More details are provided in the supporting material.
We used umbrella sampling along ξ ch to compute the PMF of stalk formation in conjunction with the coarse-grained Martini lipid force field. Figure 1D presents PMFs of stalk formation across membranes of POPC lipids using different degrees of hydration in the proximal water compartment, spanning 1.25 to 7 water beads per lipid, corresponding to 5 to 28 water molecules per lipid according to the 4:1 mapping of Martini. Here, the length of the ξ ch -defining cylinder was chosen such that ξ ch ≈ 0.2 corresponds to flat unperturbed membranes, implying that 20% of the cylinder slices are filled by lipid tail beads. ξ ch ≈ 1 corresponds to a fully formed stalk. Evidently, free energy of the stalk ∆G stalk relative to the flat membranes greatly depends on the degree of hydration, spanning values from −30 kJ/mol up to 180 kJ/mol (Fig. 1D, inset). With fewer than 2 water beads per lipid, a barrier or transition state emerges, indicating that the stalk is metastable (long-living); with less than 1.5 water beads per lipid the stalk even forms the free energy minimum. The marked dependence of ∆G stalk on the degree of hydration agrees with a recent simulation study, which used the string method together with the three-dimensional (3D) lipid tail density as order parameter to derive the minimum free energy path of stalk formation. 27 However, with the identical simulation system, the free energies for the stalk suggested by our PMFs are significantly lower (Fig. S2), possibly because our stalk state given by ξ ch ≈ 1 includes more conformational freedom than a stalk definition via the 3D density used in Ref.

27.
The Martini model allows semi-quantitative simulations, suggesting that Martini yields trends often correctly; 30  However, the free energy of the stalk with Martini 3.0beta is reduced by ∼30 kJ mol −1 relative to Martini 2.2, suggesting that the newer Martini version is more fusogenic. Simulations with four additional lipid types confirm this trend (Fig. S4). This finding suggests that the free energies reported in this study should be interpreted in terms of trends (with hydration, degree of unsaturation etc.) and not in terms of precise free energy values.

Kinetics of stalk formation agree with the PMFs
As a critical test for the quality of the reaction coordinate and, thereby, for the validity of the PMFs, we carried out free simulations of stalk formation and stalk closure (Fig. 2). We simulated a system for which ∆G stalk ≈ 0 according to the PMF (Fig. 2B). Among a total simulation time of 800 µs, we observed 8 transitions of stalk opening 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. The excellent agreement with the PMF suggests that the PMF is not affected by hysteresis problems but reports the correct ∆G stalk as given by the Martini force field. Using transition state theory, the rates are expected to follow    Evidently, apart from the degree of hydration, also the lipid type strongly influences the free energy of the stalk ∆G stalk and the free energy barrier for stalk formation. Specifically, saturated tails as found in DPPC (di-16:0-PC) disfavor stalk formation, whereas polyunsaturated tails in lipids such as di-18:2-PC greatly favor stalk formation (Fig. 4A). In addition, increasing the tail length may stabilize the stalk at low hydration of ≤6 water beads/nm 2 for which the stalk is energetically accessible (Fig. 4B). At increased hydration, however, for which the stalk is unstable, longer tails further destabilize the stalk. Likewise, the type of head group may influence ∆G stalk : the anionic head groups of dioleoyl-PS (DOPS) and DOPG disfavor stalk formation by ∼25 kJ/mol relative to the zwitterionic head group of DOPC (Fig. 4C). These findings cannot be explained merely the geometric shape of the lipids because PS headgroups are smaller than PC headgroups, whereas PG headgroups are larger than PC headgroups. Hence, the electrostatic repulsion between the anionic lipids plays a distinct role in destabilizing the stalk, possibly by effectively inducing a more positive spontaneous membrane curvature. The small zwitterionic PE headgroup facilitates stalk formation relative to the larger PC headgroup by up to 50 kJ/mol, as one may expect from the cone shape of PE lipids that may stabilize the strong negative curvature in the stalk. This trend is only inverted at very low hydration of ≤5 water beads/nm 2 , where stalk formation between PC membranes is more favorable than between PE membranes (Fig. 4C, yellow and magenta dots). Simulations with PE and PC membranes with alternative tails confirmed that, at most degrees of hydration, PE headgroups typically facilitate stalk formation relative to PC (Fig. S6).
Taken together, these data demonstrate that the free energy cost for stalk formation is controlled by a range of lipid properties. Namely, stalk formation is facilitated by (i) increased tail unsaturation, (ii) by longer tails at low hydration between the fusing membranes, (iii) by zwitterionic relative to anionic lipids, and (iv) by the small PE instead of the larger PC headgroup (at most degrees of hydration).

∆G stalk scales linearly with lipid concentration in lipid mixtures
Next, we investigated how the lipid type and concentration influence stalk formation in lipid mixtures. To this end, we set up addition 76 simulation systems containing POPC as reference lipid plus one type of lipid with increasing mole fraction. A constant degree of hydration with 6 water beads/nm 2 was used. In addition to the PC, PG, PS, and PE lipids considered above, and to obtain a comprehensive view on the influence of lipids on stalk free  x lipid = 0% corresponds to pure POPC. Typical PMFs are presented in Fig. S11. Evidently, mixing a second type of lipid into a POPC membrane may greatly stabilize or destabilize the stalk, while ∆G stalk depends nearly linearly on the mole fraction of the second lipid.
Here, the slope of the ∆G stalk -x lipid curves quantify the sensitivity of stalk formation upon the addition of a second lipid.
As expected, the geometric shape of the lipid, that is whether the lipid is cone-or inverted cone-shaped, affects ∆G stalk . Only small amounts of ceramide and diacylglycerol greatly stabilize stalks, compatible with their highly negative intrinsic curvature (Fig. 5A, dark blue, dark orange). On the contrary, single-chained lysoPC greatly destabilizes the stalk, compatible with its large positive curvature (Fig. 5A, orange). Single-chained fatty acids stabilize the stalk in the protonated form 16:0-COOH, however they destabilize the stalk in the deprotonated form 16:0-COO − (Fig. 5A, light blue, purple). Owing to the strongly increased pK a of fatty acids in lipid membranes relative to bulk water, the protonated, stalk-stabilizing form is predominantly present in biological membranes. 31 In agreement with the findings for single-component membranes presented above, ∆G stalk is modulated by replacing PC head groups of POPC with other head groups while maintaining the palmitoyl-oleoyl tails. Whereas replacing PC with PA has only a small effect on ∆G stalk , replacing PC with PE favors stalk formation (Fig. 5A, green, Fig. 5B, orange).
In contrast, replacing PC with the anionic PS or PG head groups disfavors stalk formation (Fig. 5B, light and dark blue). Likewise, sphingomyelin disfavors stalk formation, in particular the fully saturated DPSM. These findings are rationalized by the tendency of sphingomyelin of increasing the order of the aliphatic tails, thereby disfavoring the formation of a defect such as a stalk. In addition, the ∆G stalk values from lipid mixtures confirm that the addition of (poly-)unsaturated lipid tails facilitate stalk formation, whereas the addition of saturated tails hinder stalk formation (Fig. 5C). Finally, the addition of cholesterol strongly stabilizes the stalk, by up to 45 kJ/mol at 40 mol% cholesterol, in qualitative agreement with results from X-ray diffraction. 32 Stalk stabilization correlates with specific lipid enrichment in the stalk, except for cholesterol To shed light on the molecular mechanism by which different lipids influence ∆G stalk , we computed the lipid densities in simulations with a fully formed stalk (Figs 6A/B and S12-S15). The overall shape of the stalk is largely independent of the lipid composition, except that energetically unstable stalks are slightly thinner as compared to stalks that form the free energy minimum. However, the lipid composition strongly influences the spatial distribution of the individual lipid types. For instance, the stalk-stabilizing DLiPC is enriched at the stalk, whereas stalk-destabilizing DPSM is depleted at the stalk (Figs 6A/B, left column). To test whether these trends hold for other lipids, we quantified the enrichment of a lipid in the stalk as ρ stalk /ρ mem − 1, where ρ stalk and ρ mem denote the average lipid mass densities at the These findings suggest that cholesterol does not favor stalk formation owing to an intrinsic negative curvature of cholesterol, but more likely owing to its influence on the dehydration repulsion between the membranes.
Another outlier is diacyl-glycerol PODG, the by far most stalk-stabilizing lipid considered in this study. The lipid densities show that PODG accumulated not only inside the stalk, but mostly on top and below the stalk within the membrane core (Fig. S13). The accumulation within the membrane core is possible because the headgroup of diacyl-glycerol is only weakly polar. We hypothesize that similar effects play a role during biogenesis of lipid droplets, which contain large amounts of triacylglycerol and which are released from the ER membrane via a stalk-like intermediate. 34

Discussion
The marked difference between the outer and inner leaflet prompted us to screen free energies of stalk formation for a wide range of single-lipid membranes as well as for membranes of binary lipid mixtures. The results provided a comprehensive and quantitative fusogenicity map of lipids, as summarized in Fig. 6D. The data suggest that free energies of stalk formation are not determined by a single lipid property, but instead by a combination of molecular mechanisms: (i) Stalk formation is facilitated by increased tail disorder, as present in polyunsaturated lipids, whereas increased tail order hinders stalk formation, as common for sphingomyelin. (ii) As expected from previous studies, 35 lipids with negative intrinsic curvature promote stalk formation as they may stabilize the large negative curvature of the stalk. Such mechanism applies to protonated fatty acids, to PE lipids and ceramide owing to their small head groups, as well as to polyunsaturated lipids owing to the increased tail volume. A counter example is given by lysoPC that hinders stalk formation owing to its positive intrinsic curvature. 35  Therefore, we hypothesize that cholesterol may promote stalk formation by reducing the hydration repulsion, 38 thereby facilitating the formation of the first contacts between the fusing leaflets.
By simulating stalk formation between membranes with complex lipid composition, we found that the inner leaflet of a typical mammalian plasma membrane 4 is more fusogenic than the outer leaflet by ∼50 kJ/mol. These properties may be an adaption to evolutionary pressure: a fusogenic inner leaflet facilitates exocytosis, thereby increasing the fusion rates and reducing the required energy consumption by fusion proteins. A less fusogenic outer leaflet might hinder infection by enveloped viruses that fuse with the plasma membrane. The fusogenicity of lipids in Fig. 6D provides the molecular rationale for the distinct fusogenicities of the inner and outer leaflets. Namely, the increased fusogenicity of the inner leaflet is mainly caused by the increased number of polyunsaturated tails, mainly present as PS and PE lipids, and, to a lower degree, by the large PE content (cf. Table 1). The outer leaflet counteracts fusion mainly owing to the large sphingomyelin content, fewer polyunsaturated tails, and increased PC content relative to PE.
These conclusions were obtained with novel reaction coordinate for stalk formation, which allowed computationally highly efficient free energy calculations of stalk formation along thermodynamically reversible pathways. The PMFs converge rapidly, are not affected by hysteresis problems, and they yield the true free energy difference between the lamellar and the stalk state for the given force field, as confirmed by free simulations. One PMF required less than 7 hours on an inexpensive server equipped with a 6-core CPU and a consumer graphics card, hence allowing for high-throughput calculations. Our PMF calculations complement computationally more elaborate methods such as the string method, which is capable of finding the minimum free energy path of stalk formation without the need of identifying a good reaction coordinate. 27 The computational efficiency further relied on the use of the Martini coarse-grained force field that achieves a speedup in sampling by a factor of ∼1000 relative to atomistic simulations due to fewer particles, longer integration time step, and accelerated lipid diffusion. However, with the new reaction coordinate, free energy calculations of stalk formation with atomistic MD simulations are now within reach. This will allow us to test whether coarse-grained simulation primarily provide trends or whether they also provide quantitatively precise predictions of free energies during fusion.
While this study focused on stalk formation, it is evident that other stages of the fusion process are influenced by the lipid composition as well. Stalk widening as well as opening and expansion of the fusion pore involve highly curved membranes and are therefore controlled by lipid curvature. 10,39 Hence, future studies should quantify the role of complex lipid compositions on such later stages of fusion, in addition to lipid effects on the dehydration repulsion, which may add a free energy offset to the free energies of stalk formation computed here. 27 Further, it will be interesting to test whether lipid-protein interactions guide a local enrichment of fusogenic lipids and, thereby, further facilitate fusion.

Conclusions
Using a newly developed method for computing free energies of stalk formation ∆G stalk , we presented a comprehensive and quantitative fusogenicity map of lipids. The simulations showed that the lipid composition of membranes may modulate ∆G stalk by 100 kJ/mol or even more. A range of lipid properties facilitate stalk formation, including increased tail unsaturation, longer tails, smaller head groups, and cholesterol content. In contrast, stalk formation is hindered by anionic lipids and by lipids that increase the tail order such as sphingomyelin. We found that the lipid composition of the inner leaflet of a mammalian plasma membrane greatly favors stalk formation relative to the outer leaflet. The distinct fusogenicities of the two leaflets is likely an adaptation to physiological requirements and is mainly rationalized by the different content of polyunsaturated lipids, sphingomyelin, and PE lipids.

Methods Summary
MD simulations were carried out with Gromacs, version 2020.3, 40 and with an in-house modification of Gromacs 2018.8 that implements the harmonic restraint along the reaction coordinate ξ ch , originally introduced to study pore formation over membranes. 28,29 PMFs were computed with umbrella sampling using 19 umbrella windows and simulating each window for 200 ns. Interactions were described with the Martini coarse-grained force field version 2.2 if not stated otherwise. 41 Details on the reaction coordinate, simulation setup, parameters, and analysis are provided in the supporting material.