Nonuniversal impact of cholesterol on membranes mobility, curvature sensing and elasticity

Abstract

Biological membranes, composed mainly of phospholipids and cholesterol, play a vital role as cellular barriers. They undergo localized reshaping in response to environmental cues and protein interactions, with the energetics of deformations crucial for exerting biological functions. This study investigates the non-universal role of cholesterol on the structure and elasticity of saturated and unsaturated lipid membranes. Our study uncovers a highly cooperative relationship between thermal membrane bending and local cholesterol redistribution, with cholesterol showing a strong preference for the compressed membrane leaflet. Remarkably, in unsaturated membranes, increased cholesterol mobility enhances cooperativity, resulting in membrane softening despite membrane thickening and lipid compression caused by cholesterol. These findings elucidate the intricate interplay between thermodynamic forces and local molecular interactions that govern collective properties of membranes.

Introduction

Cholesterol as the most abundant lipid in plasma membranes with a total concentration of 40% in erythrocytes¹ plays a crucial role in the structure and function of cell membranes. Structurally, the steroid embeds into lipid membranes and by alignment with the phospholipids reduces the conformational freedom of their acyl chains². The increased acyl chain alignment results in an increased membrane thickness and concomitantly a decreased area per phospholipid, a property called the condensing effect of cholesterol^3,4. The different interaction of cholesterol with lipids differing in acyl chain length and saturation is one driving force in the formation of membrane phases characterized by a specific composition, packing density, and ordering⁵, with cholesterol being enriched within a liquid-ordered or L_o phase. The latter is frequently equated with raft domains⁶. Experimental evidence for these sterol- and sphingolipid-enriched domains is still debated, as is the sterol distribution between the leaflets of plasma membranes dramatically differing between experimental studies^7,8, with recent simulations hinting to a 28% increased cholesterol density within the extracellular leaflet⁹.

With respect to function, cholesterol is involved in the formation of signaling complexes in animal cell membranes, such as for G-protein coupled receptors or immune receptors. Examples include the dimerization of chemokine receptors with cholesterol acting as a sort of molecular glue¹⁰, or the formation of membrane domains, e.g. membrane rafts, that accompany the formation of the immunological synapse¹¹ and were reported to contribute to the localization of Fcγ receptors and thus immune signaling^12,13.

How does cholesterol affect the collective properties of membranes? Both in experiments^14,15,16 and simulations^17,18,19, an increase of membrane stiffness was noted upon the addition of cholesterol, in particular for bilayers containing saturated or mono-unsaturated lipids¹⁶. Differently, the bending elasticity for di-unsaturated 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) lipids was reported to be hardly affected or even decreased by the addition of cholesterol even up to a concentration of 50%, using electrodeformation analysis of giant vesicles²⁰ or X-ray scattering^16,21. Also for sphingomyelin below the main phase transition, the bending elasticity was shown to be decreased for increased cholesterol content²⁰.

At variance with this nonuniversal effect of cholesterol on the bending modulus of membranes differing in acyl chain saturation¹⁶, and combining results from neutron spin-echo experiments (NSE), solid-state 2H-NMR spectroscopy (ssNMR), and molecular dynamics (MD) simulations, Chakraborty et al. reported an increase of membrane bending rigidity with increasing cholesterol concentration also for DOPC, similar to the case for saturated membranes¹⁹. At 50% cholesterol concentration, the bending modulus showed an even 3-fold increase compared to pure DOPC membranes.

Both ssNMR and NSE measure relaxation rates^22,23. This temporal decay of observables was linked to time-averaged observables, here the membrane bending modulus, by making specific assumptions on membrane characteristics, including the approximation of a homogeneous membrane composition. The discrepancy to previous experiments addressing membrane bending fluctuations was explained in terms of a hierachical energy landscape, suggested to result in a cholesterol-dependent increased bending rigidity only on short time and length scales¹⁹. I.e., at short wave length or short relaxation times, the slow lipid diffusion entails an increase in the membrane bending modulus. However, the definition of the Helfrich bending modulus does not extend to very short length scales²⁴. It is also unclear how a cholesterol-induced membrane rigidification at short length scales may be reconciled with a membrane softening effect on larger length scales.

Interestingly, and in apparent contradiction to the increased stiffness of cholesterol-containing membranes or domains, bending of cholesterol-rich ordered domains (L_o-domain) is preferred as compared to bending of disordered, cholesterol-free domains (L_d-domain) with lower bending elasticity^25,26,27. Cholesterol flip-flops were suggested to contribute to the deformation by relaxing tension^25,27,28,29. Indeed, simulations of lipid sorting in plasma membrane tethers showed an increased concentration of cholesterol in negatively curved membrane leaflets³⁰ that follows its negative spontaneous curvature^9,31. These findings rather hint towards a decreased (apparent) bending modulus at short length scales, in disagreement to the reported increase deduced from ssNMR, NSE, and MD simulations. Interestingly, a recent atomistic MD study³² suggested a lipid-specific effect of cholesterol on bending rigidity, in agreement with electrodeformation, X-ray scattering, and tube-aspiration measurements¹⁵.

In summary, although the influence of cholesterol on the structure of lipid membranes is well investigated and described on a molecular basis, its effect on membrane elasticity and the dependency on the degree of lipid saturation is hardly understood. Given the substantial ratio of di- and poly-unsaturated lipids and of cholesterol in plasma membranes¹, the effect of cholesterol for (local) membrane rigidity is, however, of huge biological relevance for numerous processes at the plasma membrane that necessitate membrane remodeling.

Here, we demonstrate that the impact of cholesterol on membrane properties is far from universal. Molecular dynamics simulations unveil a complex interplay between cholesterol distribution, membrane curvature, and elasticity. While cholesterol consistently induces membrane thickening and lipid condensation across various lipid membranes, it exerts contrasting effects on membrane stiffness. Specifically, fully saturated lipid membranes experience significant stiffening in the presence of cholesterol, whereas membranes composed of unsaturated lipids, such as DOPC, display only weak stiffening or even softening. This softening phenomenon can be attributed to cholesterol’s negative spontaneous curvature and its increased mobility within unsaturated membranes, both laterally and transversally. These insights have profound implications for biological processes, as reduced bending elasticity facilitates localized membrane remodeling in processes like budding, endocytosis, exocytosis, and cellular signaling.

Results

Cholesterol enhances membrane fluctuations

The effect of cholesterol on lipid bilayers differing in headgroup chemistry and acyl chain saturation (1,2-dipalmitoyl-sn-glycero-3-phosphocholine, DPPC; N-(hexadecanoyl)-hexadecasphing-4-enine-1-phosphocholine, DPSM; 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine, POPC; DOPC) was analyzed for lipid bilayer (short e.g. as i:DPPC) and lipid bicelle configurations (see Fig. 1a, b for DOPC bicelles, short d:DOPC) based on molecular dynamics simulations at coarse-grained (open symbols in Fig. 1) and additionally all-atom resolution (for DPPC, DOPC, and POPC, filled symbols; see Supplementary Table 1 for a list of studied systems). Similar to previous both computational and experimental studies, a condensing effect of cholesterol was observed for all studied lipid membranes, i.e. the area per lipid decreased with increasing cholesterol concentration (Fig. 1c). The increased lipid density and the packing of cholesterol molecules between the phospholipids results in an increased order of the acyl chains (exemplarily shown for all-atom DOPC model in Fig. 1b) and an increased thickness of the lipid membranes (Fig. 1d).

Fig. 1: Lipid bilayer characteristics in the presence and absence of cholesterol.

Panel (a) presents snapshots of all-atom (AA) DOPC bicelles, without and with cholesterol molecules. DOPC within the central bicelle domain are colored in yellow (brown in bicelle rim domain), cholesterol in green, and DTPC molecules in blue. Rendering of the images was done with Blender v3.4.1⁸⁴ and Molecular Nodes v2.4.1⁸⁵. The study analyzes various features of lipid bilayers, including acyl chain order (b), area per lipid (c), and membrane thickness (d), using different setups (CG bicelles as circles and infinite, periodic bilayers as squares) and resolutions (open symbols for CG resolution and filled symbols for all-atom resolution). Error bars were smaller than the symbol size and therefore excluded.

In addition to changes in structure, addition of cholesterol affects as well the collective properties of lipids in membranes. Here, we analyzed the effect of cholesterol on membrane thermal fluctuations in the presence and in absence of cholesterol in terms of the mean membrane curvature $\bar{H}(t)$ within a central circular disc domain (radius 3 nm, see Supplementary Methods).

Fig. 2 shows the progression and distribution of mean membrane curvatures exemplarily for symmetric monovalent fully-saturated DPPC (a,b) and doubly-unsaturated DOPC systems (c,d) in the presence (blue lines) and in the absence of cholesterol (black lines, simulation lengths 4 µs). The curvature distributions could be well fitted by Gaussian distributions centered around a vanishing total curvature, differing, however, in the distribution widths (standard deviation σ in Fig. 2b, d). The results obtained did not differ significantly between infinite bilayers and bicelle systems (see Supplementary Figs. 4–8).

Fig. 2: Coupling between membrane curvature and cholesterol distribution.

Mean curvature values within circular domains of radius 3 nm were analyzed in absence (black) and presence of cholesterol (blue), separately for DPPC (at 330 K), and DOPC bilayers (320 K) at all-atom resolution. Panels a, c show the time development and panels b, d histograms for the mean curvature $\bar{H}$ with a fit assuming a Gaussian distribution. Panel e shows the low-pass filtered mean curvature values (blue) together with the asymmetric distribution of cholesterol between the two leaflets within the analyzed central circular domain (red; difference in the number of cholesterol molecules between the upper and lower leaflets, normalized to the average number within one leaflet). Panel (f) displays the time-averaged mean curvature as a function of cholesterol asymmetry for both the DPPC and DOPC bilayers.

The distribution width σ of the mean curvature values (used as a measure for the bending modulus κ_b for bicelle systems), is slightly smaller for DOPC compared to DPPC. Addition of cholesterol led to increased fluctuations for DOPC, from 0.018 nm⁻¹ to 0.021 nm⁻¹, indicative of a softened membrane. In contrast, curvature fluctuations substantially decreased for the saturated DPPC bilayer upon the addition of cholesterol, a signature of membrane stiffening.

Thus, despite a similar ordering and condensing effect of cholesterol on all investigated types of lipid bilayers, the addition of cholesterol led to enhanced fluctuations of the mean curvature for the di-unsaturated DOPC lipid bilayer, and decreased fluctuations for fully saturated DPPC membranes.

Cholesterol stiffening and softening of lipid membranes

To further quantify the effect of cholesterol on membrane fluctuations, we analyzed the membrane bending modulus κ_b for all investigated systems, i.e. bicelle and infinite systems, at CG and all-atom resolution (see Supplementary Table 4). κ_b is defined by the Helfrich Hamiltonian, which relates the membrane bending energy E_b to the local principle curvatures C₁ and C₂ by an integral over the membrane surface^33,34

$${E}_{b}=\int {dA}\left\{\frac{{{{{{{\rm{\kappa }}}}}}}_{b}}{2}{\left({C}_{1}+{C}_{2}-{C}_{0}\right)}^{2}+{{{{{{\rm{\kappa }}}}}}}_{G}{C}_{1}{C}_{2}\right\}$$

(1)

The Gaussian curvature term (κ_GC₁C₂) is disregarded in the following as it contributes a constant, the spontaneous curvature C₀ is zero for symmetric membranes. For infinite membrane systems, Fourier transformation of the lipid height field (FT-H) or, alternatively, of the lipid orientation field (FT-O)^35,36,37, and application of the equipartition theorem allows to determine bending moduli and lipid tilt moduli from fluctuation spectra³⁸ (see Methods Section for details) with high accuracy (Fig. 3c). For the bicelle, in turn, the local fluctuation method (LFM) in direct space was used⁹. It relates the fluctuation width σ in mean curvature (H = 1/2(C₁ + C₂)) averaged over the analyzed circular disc domain to the bending modulus (see Methods Section), with σ ∝ κ^−1/2.

Fig. 3: Effect of cholesterol on membrane bending modulus.

Membrane bending moduli were determined from CG MD simulations (a) and from all-atom MD simulations (b) for different system setups (bicelles as circles and infinite, periodic bilayers as squares). The bending moduli were calculated using undulation analysis (for infinite systems) and the local fluctuation method (for bicelles). Panel c shows exemplarily the spectrum for DOPC bilayer systems with 0% and 40% cholesterol for different lateral box sizes (27 nm, 55 nm, and 110 nm; $\langle|z_{{{\mathbf{q}}}} |^2\rangle$ are the amplitudes of the Fourier-transformed membrane surface). Displayed are mean values and errors as 95% confidence intervals employing parametric bootstrapping (N = 50,000 statistically independent samples) assuming Gaussian distributions of the mode-dependent amplitudes⁷⁸ (infinite bilayers) or standard errors of the mean employing block averaging (at least N = 12 independent blocks; bicelles).

At CG resolution (Fig. 3a), and in the absence of cholesterol, the bending modulus determined for infinite membrane systems is largest for DPPC (30.7 k_BT), and decreases with increasing unsaturation of the lipid acyl chains (DPSM 26.7 k_BT, POPC 22.1 k_BT, DOPC 18.8 k_BT; all CG simulations at 320 K). The simulation-derived values for κ_b are in overall good agreement to experiments (see Supplementary Table 6), FT-H and FT-O methods yielded comparable results (see Supplementary Table 4). For low cholesterol concentrations (<20 mol%), the bending moduli were mildly affected and stayed approximately constant (DPPC, DPSM) or decreased by up to 5 k_BT (DOPC, POPC). At increased cholesterol concentration of ≈40 mol%, the bending elasticity for DPPC strongly increased to >43 k_BT, i.e. the membrane substantially stiffened. In contrast, κ_b for DOPC softened by 37% compared to the cholesterol-free case, suggesting a non-universal impact of cholesterol on membrane elasticity.

The changes in bending moduli with increasing cholesterol content overall agree between bicelles and infinite bilayer systems. However, values for κ_b are consistently enlarged by 3–8 k_BT using the LFM method as compared to values obtained from undulation and orientation spectra of infinite systems.

At all-atom resolution (Fig. 3b), the stiffening effect of cholesterol on the fully saturated DPPC membrane is substantially enhanced. At a physiological concentration of 40% cholesterol, κ_b is almost fourfold increased compared to the cholesterol-free case (from 28.2 k_BT to 106.3 k_BT), in good agreement with experiments for the fully saturated DMPC^14,16. Despite a comparably high temperature of 330 K, addition of cholesterol drives the membrane into a liquid-ordered phase characterized by a significantly increased thickness (increase by 0.74 nm compared to cholesterol-free DPPC) and a very high (averaged) tail order parameter of 0.34 (Supplementary Table 3). This phase transition is insufficiently captured by the coarse-grained forcefield.

In contrast to DPPC, the bending rigidity of the di-unsaturated DOPC membrane is decreased by 31% upon the addition of 40% cholesterol at 320 K (all-atom, Fig. 3b, Supplementary Table 4), similar to the results obtained at CG resolution. However, at ambient conditions (298 K), the addition of cholesterol resulted in an increase of κ_b from 24.1 k_BT to 30.2 k_BT. This increase is comparable to results from a size distribution analysis of giant unilamellar vesicles (GUVs)³⁹.

Possible length scale-dependent effects of cholesterol-induced softening were addressed by comparing bending moduli derived from DOPC systems with lateral box lengths between 27 nm and 110 nm (Fig. 3c, CG resolution, temperature 320 K). The values of κ_b were in excellent agreement, regardless of the size of the system. That is, cholesterol induces membrane softening in DOPC bilayers even for wavelengths beyond 0.1 µm. It is important to note that usage of standard simulation parameters for neighbourlist updates, cutoffs, and restraints led to imbalances in the pressure tensor⁴⁰ and artifacts in the power spectra (not shown). These issues are particularly pronounced in simulation systems exceeding >50 nm in lateral extension, resulting, for instance, in nonphysical deformations in large coarse-grained membranes⁴⁰. In order to prevent potential biases in our spectra for both the large DOPC and DOPC/cholesterol membranes, we opted for more conservative simulation parameters, which were carefully selected and validated (Fig. 3c, Supplementary Fig. 10).

In contrast to the non-universal effect of cholesterol on membrane rigidity, an increase of the lipid tilt modulus (κ_θ) was observed for all studied lipid bilayer systems upon cholesterol addition (Supplementary Table 4), in line with the lipid condensing effect of cholesterol and the increased acyl chain ordering (Fig. 1).

Cholesterol and curvature are mutually dependent

In all studied systems, the cholesterol distribution to the membrane leaflets strongly correlated with the spontaneous curvature as exemplarily shown for the all-atom DOPC bilayer system (4 µs length) in Fig. 2e employing low-pass filtered trajectories (see Supplementary Figs 7, 8 for results of DOPC simulations at CG resolution and for bicelle systems). The appearance of temporary membrane curvature within the analyzed circular domain (radius 3 nm) implies a partial redistribution of cholesterol to the respective compressed, negatively curved leaflet. Thus, cholesterol molecules preferentially partition to regions of negative curvature. The dependency between the mean curvature and the cholesterol asymmetry is approximately linear for both all-atom and coarse-grained systems (see Fig. 2f and Supplementary Fig. 9). Between the different lipid types, it is significantly enhanced for lipids with unsaturated acyl chains (Fig. 2f). Locally, the cholesterol asymmetry may even exceed 80% for high (local) curvature (shown for i:DOPC^h:AA at 320 K, Fig. 2f and Supplementary Fig. 11). The highly dynamic coupling between cholesterol asymmetry and the local mean membrane curvature is schematically displayed in Fig. 4 for two snapshots of a DOPC all-atom bicelle: The number of cholesterol molecules between upper and lower membrane leaflet within a circular domain of 7 nm radius changed by up to 50% between the snapshots separated by 1.6 µs. This change was connected to a change in mean curvature by 0.034 nm⁻¹ between the snapshots.

Fig. 4: Curvature-dependent cholesterol dynamics and distribution.

The sketch depicts the lateral positions of cholesterol molecules within the upper and lower leaflets (upper and lower panels, respectively) of a DOPC bicelle at atomistic resolution. Snapshots at distinct times, denoted as positive curvature (a, at time t₁) and negative curvature (b, at time t₂ > t₁), showcase the dynamic redistribution of cholesterol. The central domain, highlighted in green with a radius of 7 nm, serves as a reference. Cholesterol molecules departing from the central domain of the respective leaflet between these two snapshots are depicted in red, while molecules moving into the central domain of the upper/lower leaflet are represented in orange. The inset figures were created with BioRender.com.

But how does the presence of cholesterol differentially affect membrane elasticity depending on lipid chain saturation? Cholesterol molecules are free to flip between the lipid leaflets of a bicelle and of periodic infinite membranes. The observed cholesterol flipping rates in the coarse-grained simulations were comparable for bicelles and periodic bilayers but differed between the investigated lipid types: Larger cholesterol flipping rates were observed for lipid membranes with a higher degree of acyl chain unsaturation, i.e. the cholesterol flipping rate increases in the order DPPC < POPC < DOPC, ranging between 0.4 µs⁻¹ for the fully saturated DPPC and 5.0 µs⁻¹ for the di-unsaturated DOPC (bicelle at high cholesterol concentration, see Supplementary Table 5). In addition to the increased rate of cholesterol flip-flops, the higher degree of acyl chain unsaturation as well results in increased lateral diffusion coefficients for cholesterol (see Supplementary Tables 2 and 3). The lateral cholesterol diffusion coefficient increases from (0.22 ± 0.01) × 10⁻⁶ cm²s⁻¹   for DPPC to (0.48 ± 0.01) × 10⁻⁶ cm²s⁻¹   for DOPC (CG resolution, at high cholesterol concentration), and, at all-atom resolution, from (0.09 ± 0.01) × 10⁻⁶ cm²s⁻¹  for DPPC (330 K) to (0.15 ± 0.01) × 10⁻⁶ cm²s⁻¹ for DOPC (320 K).

Our results suggest that the subtle alterations in membrane elasticity observed in certain experiments with di-unsaturated DOPC following cholesterol addition^{16,20,41,42,43}, in contrast to experiments and simulations on mono-unsaturated or fully saturated membranes (see Supplementary Table 6), can be attributed to an accelerated flipping of cholesterol molecules between leaflets and an increased cholesterol lateral diffusion coefficient within the leaflets. These dynamics are oriented towards elevating cholesterol concentration within negatively curved domains of lipid membranes. This dynamic redistribution of cholesterol, favoring locally enhanced curvatures, may potentially outweigh the anticipated stiffening effects resulting from cholesterol-induced lipid condensing and membrane thickening, as evidenced in the case of DOPC at an elevated temperature of 320 K.

This finding is additionally supported by results from MD simulations with artificially restricted mobility of cholesterol molecules: Blocked cholesterol flip-flops and restricted lateral diffusion of cholesterol led to a substantial increase of the bending moduli for both DPPC and DOPC lipid membranes (see Fig. 5 and Methods Section for methodological details). When the dynamics of cholesterol was restricted, the DOPC bending modulus increased from κ_b = 12.0 k_BT (free cholesterol dynamics) to κ_b^noflip = 15.9 k_BT (inhibited flipping of cholesterol) and κ_b^restr. = 26.4 k_BT (no flipping & restrained lateral diffusion). Similarly, the membrane stiffness of cholesterol-containing DPPC membranes was substantially increased for restrained cholesterol dynamics (from 43.0 k_BT to 71.6 k_BT). These results unambiguously show that both membranes, saturated and un-saturated, are stiffened by the structural integration of cholesterol and substantially softened by cholesterol dynamics.

Fig. 5: Effect of cholesterol dynamics on membrane bending modulus.

The bending modulus was analyzed for DPPC (a) and DOPC bilayers (b) based on unperturbed coarse-grained MD simulations in the absence of cholesterol (gray bars), in presence of 40% cholesterol, with blocked cholesterol flipping between the layers (noflip), and for simulations with blocked cholesterol flipping and restrained lateral cholesterol diffusion (restr.). Displayed are mean values and errors as 95% confidence intervals employing parametric bootstrapping (N = 50,000 statistically independent samples) assuming Gaussian distributions of the mode-dependent amplitudes⁷⁸.

Discussion

Here, we addressed the intricate interplay of cholesterol asymmetry and membrane curvature, elucidating the impact of dynamic cholesterol redistribution on the membrane bending modulus κ_b. Our simulations uncovered apparent contrasting effect of cholesterol on membrane properties. Although cholesterol consistently induces membrane thickening and lipid condensing irrespective of lipid acyl chain saturation, cholesterol exerts both a membrane stiffening and softening: A substantial stiffening was observed for fully saturated DPPC and DPSM membranes. In contrast, the stiffening was more moderate in mono-unsaturated POPC and di-unsaturated DOPC membranes at ambient temperature (all-atom resolution). Intriguingly, this effect reversed at an elevated temperature of 320 K, at both coarse-grained and all-atom resolution, revealing a significant cholesterol-induced softening of DOPC membranes.

The (local) cholesterol distribution between the membrane leaflets was observed to be strongly coupled to fluctuations of the local membrane curvature, closely mirroring the local curvature on a sub-microsecond to microsecond timescale (Fig. 2e, f, Supplementary Fig. 11). The preference of cholesterol for membrane regions with negative curvature is in line with previous observations in experiments⁴⁴ and simulations^9,30,31. Notably, the coupling strength between local membrane curvature and cholesterol redistribution was markedly heightened in di-unsaturated DOPC membranes compared to fully saturated DPPC membranes (Fig. 2f). This heightened coupling correlates with cholesterol mobility, encompassing both lateral and transversal movements through flip-flops, the latter revealing a nearly tenfold increase in flipping in DOPC compared to DPPC membranes. Interestingly, recent research on giant unilamellar vesicles suggested that differences in lipid spontaneous curvature connected to lateral lipid diffusion could lead to membrane softening⁴⁵. Our results show that a similar diffusional softening⁴⁵, here by cholesterol, may even outweigh the effects of cholesterol-induced condensation and membrane thickening on membrane elasticity.

Our findings on the effect of cholesterol on the membrane bending modulus align well with earlier experimental studies directly assessing membrane elasticity. These studies, including lipids sorting to nanotubes⁴¹, fluctuation analysis, vesicle electrodeformation²⁰, and X-ray scattering on bilayer stacks²¹, consistently reported either unchanged or slightly decreased or increased bending elasticity of DOPC upon cholesterol addition (see Supplementary Table 6 for a comparative overview). However, recent ssNMR and NSE experiments, along with MD simulations, have introduced a contrasting and controversially discussed viewpoint by suggesting a notable up to threefold increase in the DOPC bending rigidity following cholesterol addition¹⁹. Here, we demonstrate that increased lipid order parameters do not necessarily indicate an increase in membrane stiffness, as assumed in ssNMR analysis based on the assumption of a homogeneous membrane model. Our simulations show an increase in lipid order parameters but a simultaneous only weak increase (at 298 K) or even decrease (at 320 K) in the bending modulus when cholesterol is added to DOPC.

Noteworthy, similar to ssNMR and NSE, a 2–3 fold increase in DOPC bending moduli upon the addition of cholesterol was as well derived from all-atom MD simulations¹⁹. However, the real space fluctuation (RSF) analysis heavily overestimates the bending moduli for cholesterol-rich membranes (see Supplementary Table 4). RSF deduces the bending modulus from the local splay of lipid pairs separated by no more than 1 nm. The method was shown to produce results that were in reasonable agreement with Fourier-based methods and experiments for single-component membranes, with the largest deviation observed for the doubly-unsaturated DOPC (35%)¹⁷. The increase in deviation from Fourier-based approaches with increasing cholesterol concentration suggests that analysis of local structure alone may not grasp the effect of diffusional softening. Differently, and consistent with our results, a recent all-atom MD study³² showed a moderate increase of bending moduli for POPC and DOPC membranes at 300 K upon cholesterol addition (30%), by application of external forces to enforce a sinusoidal shape across periodic membrane systems combined with an umbrella sampling protocol⁴⁶. However, the study concludes that (i) membrane stiffening is not affected by cholesterol on length scales below 10 nm and (ii) did not observe and, therefore, excludes local cholesterol enrichment as a cause of membrane stiffening. Our analysis of local mean curvature distributions for circular regions (radius 3 nm) unambiguously reveals that cholesterol addition does have an effect on curvature generation and membrane rigidity also on small length scales (Fig. 2a-d), and that curvature generation and leaflet cholesterol asymmetry are tightly interconnected (up to 80% cholesterol asymmetry in circular areas of radius 3 nm, Fig. 2f; Supplementary Fig. 11).

The simulation data implies that cholesterol adopts multiple roles in the mechanics of cellular membranes: Its lipid condensing effect stabilizes plasma membranes with respect to permeability⁴⁷ or membrane rupture⁴⁸. This effect is pronounced when cholesterol induces a phase transition toward the liquid-ordered phase. In turn, the mobility of cholesterol in unsaturated membranes together with its spontaneous negative curvature counteracts the stiffening induced by lipid condensation and maintains a low or even reduced bending rigidity compared to that of cholesterol-free membranes. Yet, the effect of cholesterol on the mechanical properties of plasma membranes remains enigmatic, given its asymmetric lipid distribution with twofold more unsaturated and more mobile lipids in the cytoplasmic leaflet, as reported for red blood cell (RBC) membranes^1,49, and likely heightened cholesterol density in the exoplasmic leaflet⁹. Intriguingly, recent studies indicate a remarkably low bending modulus for cholesterol-rich RBC membranes of only 4–6 k_BT⁵⁰, contrasting with an anomalous stiffening reported for cholesterol-free model membranes with asymmetric lipid compositions^51,52. Taken together, cholesterol likely contributes to reconciling the seemingly conflicting roles of plasma membranes, i.e., creating an effective barrier against small molecules and ions while concurrently ensuring the elasticity and plasticity essential for processes requiring dynamic membrane remodeling.

Our findings propose that the induction of local curvature by transmembrane proteins or membrane-associated proteins is facilitated by the local redistribution of cholesterol, effectively alleviating membrane stress. Notably, we observed an asymmetric distribution of cholesterol between membrane leaflets in the proximity of a curvature-inducing voltage-sensing potassium channel⁹. This effect is possibly amplified by the reported enrichment of polyunsaturated lipids near membrane proteins^53,54 that allow enhanced cholesterol diffusion and flip-flop rates.

Intriguingly, cholesterol depletion was found to shift the activation pressure of the curvature-sensing PIEZO1 channel to higher values⁵⁵. Given that PIEZO1 activation involves substantial local membrane deformation, transitioning from a highly curved conformation to a flattened structure⁵⁶, cholesterol depletion and the resultant possibly increased membrane rigidity may heighten the energetic cost of such transitions. Furthermore, the negative spontaneous curvature of cholesterol, along with its impact on membrane rigidity, is anticipated to play a pivotal role in membrane fusion and fission events, i.e. during endo- and exocytosis. For instance, cholesterol has been identified as a pathway switch for influenza virus infection of cells by lowering the energy required for forming a hemifusion stalk⁴⁸. Similarly, cholesterol depletion has been reported to impede curvature generation in clathrin-mediated endocytosis⁵⁷. These insights underscore the intricate interplay between cholesterol, local membrane curvature and membrane rigidity, and membrane remodeling processes, shedding light on the impact of local membrane heterogeneity related to cholesterol mobility and its negative spontaneous curvature on collective membrane properties in cellular events.

In summary, our simulation study reflects similar to earlier experiments and simulations a condensing effect of cholesterol on membranes that is independent of the particular lipid composition. Differently, cholesterol has a non-universal effect on the stiffness of membranes: It leads to a substantial increase of bending rigidity for saturated lipid membranes, but only weak stiffening or even softening for membranes composed of unsaturated lipids such as DOPC. The observed softening could be traced back to the increased cholesterol mobility in unsaturated membranes, both lateral and transversal. The cholesterol-induced changed membrane mechanical properties have far-reaching implications for biological processes, the reduced bending elasticity will locally ease membrane remodeling processes involved e.g. in budding, endo- or exocytosis, and signaling.

Methods

Molecular dynamics simulations

Several different lipid types were studied at different cholesterol concentrations in a lipid bicelle setup (denoted as e.g. d:DPPC) and in infinite membrane systems (denoted as e.g. i:DPPC). The chosen bicelle setup makes use of short-chained 1,2-ditridecanoyl-sn-glycero-3-phosphocholine (DTPC) lipids and flat-bottomed potentials to reduce the line tension, vesiculation and free diffusion across the bicelle rim, issues usually associated with bicelle systems⁹. Importantly, the bicelle model is characterized by largely unbiased thermal fluctuations as opposed to standard infinite bilayer systems⁹ (see Supplementary Figs 2–8 for comparisons of mean curvatures and cholesterol asymmetries for bicelles and infinite bilayers). All systems were simulated using the GROMACS 2020(2021, 2023) simulation package⁵⁸ at coarse-grained (CG) resolution employing the MARTINI force field (version 2.2)^59,60,61 (see Supplementary Fig. 1 for studied lipids in MARTINI representation), and, for selected systems, additionally at all-atom (AA) resolution using the CHARMM36 force field⁶². Details of the system setup are provided below and in Pöhnl et al.⁹. Details on the different lipid compositions and simulation lengths for bicelle sytems and infinite lipid bilayer systems are provided in Supplementary Table 1.

Setup of simulation system and simulation parameters

Bicelle systems

The bicelle systems consist of a rim region and a central region with defined lipid composition to investigate membrane characteristics. First, a circular membrane patch of ≈9.2 nm was setup based on an initial hexagonal lipid bilayer constructed with insane⁶³, subsequent solvation in a larger rectanglular box and simulation for 20 ns with position restraints (force constant of 1,000 kJ/mol/nm²) in z-direction on the phosphate beads. Next, two rings of rim lipids (≈380 DOPC; ≈370 DTPC) were placed around the circular membrane patch, resolvated (≈120,000 CG water molecules), and equilibrated for 20 ns with position restraints on lipids of the circular membrane patch. For production simulations, cylindrical flat-bottomed potentials were applied on the phosphate beads to restrict the motion of lipids within the central region to a circular domain of radius 10.5 nm, and to restrain rim lipids to the region beyond a radius of 8.1 nm (force constant 500 kJ/mol/nm²). All systems were studied at CG resolution using the MARTINI force field^59,60,61 with in total ≈135,000 CG beads and a lateral box size of ≈ 32 nm. The temperature was kept at 320 K via the v-rescale algorithm⁶⁴, and the Parrinello-Rahman algorithm^65,66 with a compressibility of 3·10⁻⁴ bar⁻¹ was used for semi-isotropic pressure coupling (isotropic for bicelle systems). A reaction field algorithm with a cut-off of 1.1 nm and a relative dielectric constant of 15 was used for electrostatic interactions. Lennard-Jones interactions were treated with a single cut-off of 1.1 nm and the time step was 20 fs.

For selected systems, all-atom simulations were performed using the CHARMM36⁶² force field. Initial structures were obtained by transferring equilibrated and approximately flat bicelle structures to atomistic resolution using the initram protocol⁶⁷. Resulting systems had ≈1 Mio atoms. For bicelle systems, flat-bottomed potentials were applied on the phosphate atoms of the obtained stuctures using the same radii and force constants as in the coarse-grained systems. The temperature was kept at 320 K via the v-rescale algorithm⁶⁴. Isotropic pressure coupling to 1 bar was applied using the Parrinello-Rahman algorithm^65,66 with a compressibility of 4.5 · 10⁻⁵ bar⁻¹. Van der Waals forces were smoothly switched to zero (between 1.0 nm and 1.2 nm) and the electrostatic interactions were treated with the PME method⁶⁸. The time step was 2 fs. The all-atom systems were studied at a salt concentration of 0.15 M NaCl.

Infinite membrane systems

Infinite lipid bilayer systems at CG resolution were created with insane⁶³, further detailed in⁹. The typical lateral box size was ≈27 nm, the membranes were solvated with ≈75,000 CG water beads, resulting in systems with ≈100,000 CG beads. DOPC membranes at 0% and 40% cholesterol were additionally studied for four-fold and 16-fold increased lateral sizes, the box sizes were ≈55 nm and ≈110 nm, respectively. In order to avoid possible systematic imbalances in the pressure tensor due to missed long-range attractive van der Waals interactions as reported recently [40], we chose different settings for neighbour list updates (update every 10 steps instead of 20), an increased outer cutoff radius of r_l = 1.5 nm, as well as modified settings for the LINCS algorithm to solve bonded constraints (lincs_order = 12, lincs_iter = 2). These constraint settings were shown to avoid temperature gradients for cholesterol-containing membranes⁶⁹. We observed that both the standard neighbour list settings as well as the default settings for LINCS (lincs_order = 4, lincs_iter = 1) result in artificial membrane deformations for very large systems (>50 nm) that are reflected in artificially enhanced undulations for low q values / large wavelengths (not shown).

Infinite all-atom DOPC (DPPC) systems at 0% and 40% cholesterol at a temperature of 320 K (330 K) were setup using either backmapping or the CHARMM-GUI^70,71 (lateral box sizes between 25 nm and 28 nm, ≈900,000 atoms, ≈170,000 water molecules). DOPC membranes were additionally studied at temperatures of 310 K and 298 K (0% and 40% cholesterol, lateral box size ≈17 nm, ≈300,000 atoms, ≈60,000 water molecules).

Restrained cholesterol dynamics

The role of cholesterol mobility and distribution for membrane elasticity was additionally addressed for infinite lipid bilayers (lateral box size ≈18 nm, ≈33,000 CG beads, ≈20,000 water beads) at CG resolution with blocked cholesterol flipping and restricted lateral cholesterol diffusion. Starting from a homogeneous cholesterol distribution of 40 mol% in either DOPC or DPPC lipid bilayers, the cholesterol motion was restricted in lateral direction around the initial position by cylindrical flat-bottomed potentials with a radius of 0.5 nm (force constant of 5,000 kJ/mol/nm²). Cholesterol flips were prevented by adding a short-ranged repulsive interaction between the ROH bead in cholesterol and the last beads of the lipid fatty acyl chains (Lennard-Jones parameters c₆ = 0 kJ mol⁻¹nm⁶; c₁₂ = 1 kJ mol⁻¹nm⁶).

Analysis

Custom analysis code was written using the MDAnalysis library^72,73. Membrane properties of lipid-bicelle systems were analyzed on the central bicelle domain (radius 7 nm). For infinite systems, the whole system was considered. Unless otherwise specified, errors were estimated using block averaging on the equilibrated part of the trajectory. The number of independent blocks was determined by monitoring the convergence of the standard error. Both for coarse-grained and all-atom simulations at least several microseconds were used for the analysis of all properties. Equilibration times were chosen such that any drift in the area per lipid, thickness, and order parameters vanished. The membrane properties are summarized in Fig. 1, and in Supplementary Tables 2 and 3.

Membrane surface

The PO4 beads (P atoms for atomistic systems) were used to calculate the surfaces of the membrane leaflets. Surfaces of the upper z^u(x, y) and lower z^l(x, y) leaflets were defined as height values on a 2D lateral grid (see also^74,75). Grid cells (i, j) get assigned the weighted sum of the z-positions of surrounding PO4 beads (see Supplementary Fig. 1) as z_ij-value, using Gaussian weights $w=\exp \left(-{d}_{{ij}}^{2}/2{{{{{{\rm{\sigma }}}}}}}^{2}\right)$ with d_ij as the lateral distance between the grid point (i, j) and the PO4 bead. The standard deviation was chosen to σ = 0.8 nm and the grid spacing to ∆x = ∆y = 0.5 nm. The bilayer surface was defined as the mean of the monolayer surfaces z_ij = (z^u + z^l)/2 and the local membrane normals by ${N}_{{ij}}^{u/l}=\pm \left(-{\partial }_{x}{z}_{{ij}}^{u/l},-{\partial }_{y}{z}_{{ij}}^{u/l},1\right)$.

Local membrane curvature

The bilayer surface z_ij was used to calculate the local membrane curvature:

$${H}_{{ij}}=\frac{\left(1+{\left({\partial }_{x}{z}_{{ij}}\right)}^{2}\right){\partial }_{y}^{2}{z}_{{ij}}-2{\partial }_{x}{z}_{{ij}}{\partial }_{y}{z}_{{ij}}{\partial }_{x}{\partial }_{y}{z}_{{ij}}+\left(1+{\left({\partial }_{y}{z}_{{ij}}\right)}^{2}\right){\partial }_{x}^{2}{z}_{{ij}}}{{2\sqrt{1+{\left({\partial }_{x}{z}_{{ij}}\right)}^{2}+{\left({\partial }_{y}{z}_{{ij}}\right)}^{2}}}^{3}}$$

(2)

Local curvature values H_ij(t) were averaged ($\bar{H}\left(t\right)$) over all points i, j within a circular region around the bicelle center with radius R. For infinite systems, a similar circular region was used to calculate $\bar{H}\left(t\right)$.

Order parameter

For coarse-grained systems the order parameter

$${P}_{2}=\frac{1}{2}\left\langle 3{\cos }^{2}\left({{{{{\rm{\theta }}}}}}\right)-1\right\rangle$$

(3)

was calculated for bond vectors in the lipid tail. Here θ is the angle between the bond vector and the local membrane normal in the upper, respectively lower, leaflet ${N}_{{ij}}^{u/l}$. The local normals were assigned from the laterally closest grid cell to the center of mass of the C1A and C1B beads (MARTINI nomenclature, see Supplementary Fig. 1). The expectation value was taken over all lipids and all tail bonds.

For all-atom systems, the order parameter $-{S}_{{CD}}=\frac{1}{2}\left\langle 3{\cos }^{2}\left({{{{{\rm{\theta }}}}}}\right)-1\right\rangle$ was calculated for all carbon atoms of each chain. Here θ is the angle between the local membrane normal and the C-H bond vectors. The local normal was calculated based on all lipids within a distance of 1.8 nm from the respective lipid. Order parameters for all investigated systems are provided in Supplementary Tables 2 and 3.

Membrane thickness

The thickness was calculated for each lipid position within the central analysis domain and separately for both membrane leaflets. PO4 MARTINI beads (P atoms for AA simulations) were taken as lipids positions. First, the local normal was determined based on lipids within 1.8 nm distance of the analysis lipid. The (local) thickness was calculated as the distance along the local membrane normal between the analysis lipid and the averaged (projected) positions of lipids of the opposing leaflet within a distance of 5.5 nm of the analysis lipid. Values provided in Supplementary Tables 2 and 3 are averages over all analyzed frames and all lipids.

Area per lipid

Local Voronoi tesselation was performed for each lipid (see also⁷⁶). The lipid positions were defined as the centers of mass of the GL1 and GL2 beads (AM1 + AM2 for DPSM; ROH for CHOL; glycerol oxygen atoms for lipids and hydroxyl oxygen atom for cholesterol in AA systems). For each lipid, lipid positions closer than 4 nm were considered. The area per lipid was calculated as the area of the central Voronoi cell perpendicular to the local membrane normal and averaged over all analyzed frames and lipids (see Supplementary Tables 2 and 3).

Lipid diffusion

The Einstein relation in two dimensions

$$D=\mathop{{{{{\mathrm{lim}}}}}}\limits_{t\to \infty }\frac{MSD(t)}{4t}$$

(4)

relates the lateral mean square displacement MSD(t) to the diffusion coefficient D. Lateral displacements were calculated in the xy-plane for the PO4 beads (P atoms for AA systems; ROH, respectively hydroxyl oxygen, for cholesterol) and the diffusion coefficient D fitted for lag times between 5 ns and 10 ns. Analysis was restricted to lipids initially within a central circular analysis domain of radius 3 nm.

Leaflet assignment and flipping rate of cholesterol

For every frame, all cholesterol molecules were assigned to the closest membrane leaflet (distance between the ROH bead, hydroxyl oxygen for AA systems, and PO4 beads, P atoms for AA systems). The cholesterol assignment was also used to calculate flipping rates. Flips were counted only if the cholesterol molecule stayed for ≥10 ns within the new leaflet.

Membrane elasticity

The membrane bending modulus κ_b was assessed employing four different methods: the spectrum of membrane undulations (Undulation Spectrum Method)^35,38,77,78, the spectrum of lipid orientations (Orientation Spectrum Method)^35,38,77,78, using the probability distribution of local membrane curvatures (Local Fluctuation Method)⁹, and the distribution of splay values of neighbouring lipid pairs (Real Space Fluctuations Method)^79,80.

The Undulation Spectrum Method uses the spectrum of the Fourier-transformed membrane surface as defined by the lipid head positions (z_q) to calculate the bending modulus κ^u, the effective tilt modulus κ_θ, and the soft-mode divergence q_c³⁸:

$$\left\langle|z_{{{\mbox{q}}}} |^2 \right\rangle=\frac{1}{1-(q/q_c)^2} \left( \frac{ k_{{{\mbox{B}}}}T}{\kappa_b^{{{\mbox{u}}}} q^4}+\frac{ k_{{{\mbox{B}}}}T}{\kappa_\theta q^2} \right)$$

(5)

For CG systems, the center of mass of the PO4, GL1, and GL2 beads (PO4, AM1, and AM2 for DPSM) was used for the lipid head position. For AA systems, the lipid head was defined as the center of mass of the headgroup phosphorus atom (P) and the glycerol backbone carbon atom (C2).

The Orientation Spectrum Method uses the spectrum of the lipid director field in Fourier space (n(q)) to obtain the bending modulus ${\kappa }_{b}^{{{{{{\rm{o}}}}}}}$:

$$\left\langle|{{\mbox{n}}}_{{{\mbox{q}}}}^{\parallel} |^2 \right\rangle=\frac{1}{1-(q/q_c)^2} \frac{ k_{{{\mbox{B}}}}T}{\kappa_b^{{{\mbox{o}}}} q^2}$$

(6)

Lipid directors are defined as vectors connecting lipid head and lipid tail. Heads are defined as above. For CG systems, the center of mass of the last beads of the fatty acyl chains was used for the lipid tail position. For AA systems, the tail was defined as the center of mass between the terminal methyl carbons of the fatty acyl chains.

The 95% confidence intervals for the moduli obtained through Fourier methods were obtained using parametric bootstrapping following the approach by by Ergüder et al.⁷⁸. To estimate the standard errors (δ|z_q | ²) for the spectrum amplitudes (|z_q | ²), block averaging was used. Subsequently, 50,000 new spectra were generated by sampling amplitudes from a Gaussian distribution centered at |z_q | ² with a standard deviation of δ|z_q | ². The confidence intervals were determined by fitting moduli to these generated spectra.

Since Fourier Methods are not applicable for finite bicelle systems, the Local Fluctuation Method (LFM) was employed instead. LFM directly relates local fluctuations in membrane curvature H to the bending modulus κ_b. Briefly, it explicitly evaluates the membrane-bending energy

$${E}_{b}=\int {dA}2{{{{{{\rm{\kappa }}}}}}}_{b}{\left(H-{H}_{o}\right)}^{2}$$

(7)

within a defined circular area of the membrane on a grid, assuming a vanishing spontaneous membrane curvature H_o, and local deviations ${\left(\Delta {H}_{i}\right)}^{2}={\left({H}_{i}-\bar{H}\right)}^{2}$ and the degeneracy of states to be independent of the mean curvature $\bar{H}$.

The probabilities of mean curvatures $p\left({\bar{H}}_{1}\right)$ and $p\left({\bar{H}}_{2}\right)$ are then related by

$$\frac{p\left({\bar{H}}_{1}\right)}{p\left({\bar{H}}_{2}\right)}=\exp \left(-\frac{2{{{{{{\rm{\kappa }}}}}}}_{b}A\left({\bar{H}}_{1}^{2}-{\bar{H}}_{2}^{2}\right)}{{k}_{B}T}\right)$$

(8)

Using the variance σ² of the normally distributed $p\left(\bar{H}\right)$, the bending modulus κ_b is given by:

$${{{{{{\rm{\kappa }}}}}}}_{b}=\frac{{k}_{B}T}{4A}\cdot \frac{1}{{{{{{{\rm{\sigma }}}}}}}^{2}}$$

(9)

Here no lipid tilt modulus is included, values for the bending elasticity have thus to be considered as an upper bound for the true elasticity.

For a comparative analysis of bicelle and infinite membrane systems we chose a circular analysis domain for the LFM with radius 3 nm. Circular patches of this size displayed for all systems a comparative degree of cholesterol asymmetry (see Supplementary Fig. 3) and coupling of cholesterol asymmetry to the local mean curvature (Supplementary Figs 4–8). Importantly, fluctuations of the mean curvature within the circular analysis domains showed very good agreement for infinite periodic bilayers and bicelle systems (Supplementary Figs 4–8). In addition, we compared (apparent) κ_b values derived from curvature fluctuations (LFM) to those from Fourier methods applied to the periodic bilayer systems (see Fig. 3a), revealing an overall good agreement between the different methods.

Additionally, values for the bending modulus were compared to those obtained using the Real Space Fluctuations Method (RSF)^79,80. In RSF, the distribution of splay values S = ∇n − ∇N between pairs of neighbouring lipids is used to calculate the bending modulus κ_b^R. The lipid directors n are defined as above for CG systems (for cholesterol the R1, respectively, R5 beads are used). For AA systems, the center of mass of the three terminal tail carbons of both chains are used as tails and the heads are defined as the center of mass of the headgroup phosphorous atom (P), the glycerol backbone atom (C2) and the first three tail carbon atoms (C21, C22, C23, C31, C32, C33). For cholesterol, the director connects the C3 and C17 ring-carbon atoms (atom names taken from CHARMM nomenclature).

For larger systems and hence larger undulations it was proposed before to use a local instantaneous surface instead of the averaged surface employed in the original publications⁸¹. Here, the local instantaneous normal N was calculated as described for the thickness calculation, using the center of mass of the C1A and C1B beads (R1 bead for cholesterol) for the surface definition in CG systems. For AA systems, the center of mass of the first three atoms of the tail carbons (C21, C22, C23, C31, C32, C33) was used (C3 atom for cholesterol). For all pairs of neighbouring lipids (surface distance <1.1 nm) the splay values were calculated as the covariant derivative of n − N orthogonal to the local membrane normal and in the direction that connects the surface atoms. For each pair of lipid types i and j the distribution of splay values is described by the Boltzmann distribution P(S):

$$P(S)=C \exp\left(- \frac{\chi_{ij} S^2 A_l}{2 k_{{\mbox{B}}}T } \right)$$

(10)

The splay moduli χ_ij for lipid species pairs ij are obtained from a quadratic fit to −2k_BT ln (P (S)) /A_l, where A_l is the area per lipid. Monolayer bending moduli κ_m^R are calculated from the splay moduli χ_ij:

$$\frac{1}{\kappa_m^{{{\mbox{R}}}}}=\frac{1}{\phi} \mathop{\sum}\limits_{i,j} \frac{\phi_{ij}}{\chi_{ij}}$$

ϕ_ij is the number of neighbour pairs of type ij and ${{{{{\rm{\phi }}}}}}={\sum }_{i,j}{{{{{{\rm{\phi }}}}}}}_{{ij}}$. The bilayer bending modulus ${{{{{{\rm{\kappa }}}}}}}_{b}^{{{R}}}$ is the sum of the monolayer bending moduli ${{{{{{\rm{\kappa }}}}}}}_{m}^{{{R}}}$.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.