Molecular mechanisms and energetics of lipid droplet formation and directional budding

Abstract

The formation and budding of lipid droplets (LDs) are known to be governed by the LD size and by membrane tensions in the endoplasmic reticulum (ER) bilayer and LD-monolayers. Using coarse-grained simulations of an LD model, we first show that ER-embedded LDs of different sizes can form through a continuous transition from wide LD lenses to spherical LDs at a fixed LD size. The ER tendency to relax its bilayer modulates the transition via a subtle interplay between the ER and LD lipid densities. By calculating the energetic landscape of the LD transition, we demonstrate that this size-independent transition is regulated by the mechanical force balance of ER and LD-tensions, independent from membrane bending and line tension whose energetic contributions are negligible according to our calculations. Our findings explain experimental observation of stable LDs of various shapes. We then propose a novel mechanism for directional LD budding where the required membrane asymmetry is provided by the exchange of lipids between the LD-monolayers. Remarkably, we demonstrate that this budding process is energetically neutral. Consequently, LD budding can proceed by a modest energy input from proteins or other driving agents. We obtain equal lipid densities and membrane tensions in LD-monolayers throughout budding. Our findings indicate that unlike LD formation, LD budding by inter-monolayer lipid exchange is a tension-independent process.

---

Introduction

Lipid droplets (LDs) are cellular organelles with a specific structure, which are used by living cells to store and supply lipids for energy and lipid metabolism and membrane synthesis.^1–7 Previously perceived as passive energy packages, LDs have recently emerged as highly dynamic, mobile organelles that are actively implicated in diseases such as obesity, hepatitis, and cancer,^3,8–10 and are involved in steroid synthesis and cellular development.⁴ LDs are distinguished from typical organelles by unique features: unlike other organelles with an aqueous interior, LDs enclose an oily compartment; unlike bilayer-bounded organelles with negligible membrane tensions, LDs are surrounded by a lipid monolayer with a relatively high tension;¹¹ the initial core of LDs inside the endoplasmic reticulum (ER), i.e. the nascent LD, has a typical size of a few tens of nanometers ranging from 30 to 60 nm in diameter.^12,13 These characteristic features have significant implications for the initial formation of symmetric LDs (with respect to the ER-bilayer) and subsequent budding of asymmetric LDs from the ER.

LD formation, including initial nucleation and shape transition of symmetrical LDs, appears to be mainly driven and regulated by membrane tensions. Under excess energy or stress conditions, neutral lipids are synthesized in the inter-leaflet space of the ER membrane,¹⁴ thereby creating a tense ER-bilayer.¹⁵ Beyond a concentration threshold, these neutral lipids – initially dispersed between the leaflets of the ER-bilayer – coalesce to form nascent LDs. This demixing process,^1,16 driven by the reducing contact between the neutral and membrane lipids, helps the ER relax its tense membrane. The LDs, which comprise a hydrophobic core covered by two lipid monolayers,^13–15,17 have been experimentally observed in the form of both eye-shaped lenses (simply referred to as LD lenses hereafter) and spherical LDs,^13,16 suggesting that these structures are stable under certain conditions. Stability analysis of a model LD relates LD morphological transition (from LD lenses to spherical LDs) to the LD size, membrane bending stiffness, and an imbalance between the surface tensions experienced by the ER-bilayer and LD-monolayers.¹¹ These ER and LD-tensions, known to be related to membrane lipid density,¹⁸ thus significantly contribute to forming and maintaining ER-embedded LDs with different shapes. Nevertheless, the molecular details of how tense LDs form from initial tensionless ER-bilayers and the free energy landscape of LD formation with the particular contributions of membrane tension, membrane bending, and line tensions remain unexplained.

LD transport and degradation occur through budding of LDs from the ER toward the cytosol.^19–21 This directional budding helps LDs reach cytosolic proteins and lipid-releasing enzymes, which are crucial for LD metabolism.^12,17,20 LDs, emerged from the ER, grow in size leading to larger LDs on the micron scale, which either detach from or remain attached to the ER membrane, the latter being more likely.^17,19,20 Multiple mechanisms have been hypothesized to drive the directional budding of LDs wherein the required membrane asymmetry between cytosolic and luminal LD-monolayers is imposed by different means: membrane curvature, generated by asymmetric lipid composition^22–26 and protein insertion;^12,13 tension imbalance between the LD-monolayers;¹⁵ and lipid number asymmetry between the LD-monolayers.^12,27

The specific contributions of membrane curvature asymmetry, membrane tension imbalance, and lipid composition asymmetry to LD budding is not clear. Membrane curvature, generated in the cytosolic monolayer by inserting proteins or intrinsically-curved lipids, is shown to assist LD directional budding.^26,28 Several proteins and enzymes contribute to LD formation and its directional budding by enforcing membrane tension asymmetry, via regulating lipid composition, or like fat-storage-inducing FIT2 and ARF1 through generating asymmetric lipid number.^12,14,29,30 The particular role of other proteins such as Seipin, an integral membrane protein that forms a scaffold around the cytosolic LD-monolayer, is also not certain. While Seipin was experimentally found to redundantly take part in the LD budding,³¹ more recent experiments and coarse-grained simulations suggest that it effectively induces LD budding.³²

Regardless of the driving mechanism, LD budding toward the cytosol requires a continuous area growth in the cytosolic LD-monolayer and a corresponding area reduction in the luminal one. This idea is further supported by the fact that the LD-bounding monolayer likely originates from the ER cytosolic leaflet.^19,33 The area asymmetry can be generated in different ways such as lipid biosynthesis,^27,34 fusion of the transport vesicles to the ER,³⁴ ER lipid flippase,³⁵ conversion of neutral lipids to phospholipids by FIT2,³⁶ and lipid transfer from the luminal to the cytosolic LD-monolayer by proteins such as FIT2 and ARF1. It is not clear which of these pathways are responsible for creating membrane area asymmetry and what mechanisms are driving them.

LD budding has been explained in terms of a discontinuous second-order transition from symmetric LDs to completely-budded asymmetric LDs, regulated by the LD size, and ER and LD-tensions.^11,28 Across a wide range of LD sizes, initially stable LDs with large ER-tensions become metastable for smaller tensions and eventually unstable for negligible ratios of ER to LD-tensions, where the fully-budded LD is stable. Like formation of symmetric LDs, the budding and emergence of asymmetric LDs is thought to be governed by membrane tensions.

The small size of nascent LDs, below the standard optical microscopy resolution, and the dynamic nature of their formation, make it difficult to investigate LD formation and budding experimentally. Computer simulation of coarse-grained models with molecular resolution offers an alternative platform to dynamically monitor LD formation and budding and to elucidate mechanistic aspects of LD biogenesis by precise measurement of its mechanical properties.^37,38 Molecular dynamics simulations have been successfully used to study different aspects of LD formation such as the role of proteins like Seipin³⁹ and Caveolin,⁴⁰ implications of Triglyceride blisters for LD biogenesis,^41,42 surface properties of LDs,⁴³ and LD budding by lipid exchange and LD-tension imbalance.^12,14 Here, we used coarse-grained simulations to explore mechanistic aspects and energetics of LD shape transition and budding, which are poorly understood despite extensive experimental and modelling studies.¹

To understand LD formation, we started with identifying the relevant LD parameters by exploring spontaneous LD formation. Consistent with prior studies, membrane tensions, proportional to the lipid density or area-per-lipid, were found to determine the morphological behavior of LDs.^14,18 Contrary to earlier understanding,¹¹ we observed various LD shapes, including stable spherical LDs and LD lenses, across several LD sizes, even for very small LDs. Based on these findings, we hypothesized that stable symmetrical LDs of various shapes form via a continuous transition driven by ER and LD tensions, irrespective of LD size.

To corroborate this hypothesis and to better understand morphological transitions and stability of the symmetrical ER-embedded LDs, we then explored the transition from LD lenses to spherical LDs at fixed LD sizes. By calculating the free energies, we demonstrated that by reducing the ER-tension, a continuous LD transition takes place from early-stage LD lenses with large ER-tensions to spherical LDs with negligible ER-tensions. In line with previous studies,^15,18 a subtle interplay between the area-per-lipids in the ER and LD membranes was found to control membrane tensions and the continuous LD transition. The same interplay between the areas-per-lipids in the ER and LD membranes was also found to build large membrane tensions in LD monolayers. These findings promote the previous understanding of LD formation based on size-dependent discontinuous LD transitions.^11,28

Starting from symmetrical spherical LDs embedded in the ER, we then exchanged lipids from the luminal to the cytosolic LD-monolayer, demonstrating that the area asymmetry created by lipid exchange between the LD-monolayers can act as a novel mechanism for LD budding. Remarkably, we observed that LD budding by lipid exchange is energetically neutral due to the nearly equal areas-per-lipid and resulting equal tensions in LD-monolayers. This intriguing result suggests that LD budding can proceed in a quasi-equilibrium manner through intermediate stable states with a modest energy input. Yet this slight energy input, likely supplied by proteins, is essential to drive the LD budding. In contrast to previous models that described LD budding as a tension-induced discontinuous transition from symmetric to asymmetric LDs,^11,28 our results unveil a neutral transition from symmetric spherical LDs to fully-budded asymmetric LDs. Thus, we demonstrate that, unlike LD formation, LD budding is not driven by membrane tensions but can rather occur through lipid exchange between the LD-monolayers, driven by a modest energy input.

The shape of the nascent LDs depends not only on the ER and LD-tensions but also on the tension along the contact line of LDs with the surrounding ER-bilayer.^28,44 We demonstrated that even for the nano-sized LDs neither line tension nor membrane bending significantly contribute to the LD formation. Comparing our simulation results with the theory of membrane elasticity, we found a strong agreement between the LD shapes predicted by the theory and those observed in the simulations.

Methods

Coarse-grained simulations

We performed dissipative particle dynamics (DPD) simulations^45–47 of a coarse-grained model to study LD shape transition and budding. Our simulation method has been widely used to study various aspects of both artificial membranes⁴⁸ and cellular membranes, including membrane fusion,⁴⁹ cellular membrane rupture and damage,⁴⁵ membrane domain simulation,⁵⁰ vesicle dynamics in shear flow,⁵¹ and membrane-cholesterol interactions.⁵² The soft molecular interactions of DPD allow faster dynamics with larger time steps.^46,53,54 Therefore, DPD has been particularly effective in reproducing mechanical properties of lipid membranes^55,56 such as membrane tensions^57–59 as well as the phase behavior of membranes,^46,52 which mainly contribute to LD formation and budding.

Our coarse-grained LD model is composed of lipid membranes representing the ER-bilayer and LD-monolayers, an aqueous phase representing the luminal and cytosolic liquids, and an oil phase that represents neutral lipids i.e. the fatty core of the LDs. We note that both cytosolic and luminal liquids are taken to be identical, simply represented by water beads in our model. Coarse-grained phospholipids consist of three hydrophilic head-group beads H and two hydrophobic chains each having six tail-group beads T, see ESI,† Fig. S1a. Both cytosolic and luminal aqueous liquids are composed of water beads W. Each neutral lipid is modelled as a single chain consisting of four oil beads O, see ESI,† Fig. S1b. All simulations were performed in LAMMPS simulation package.⁶⁰

All different coarse-grained bead types H, T, W, O have a diameter of d ≈ 0.8 nm, which defines the length scale in our simulations.⁶¹ These beads interact by a conservative DPD force

where r_ij is the distance between two beads, i and j, and the unit vector r̂_ij connects the bead i to the bead j. Whereas the interaction force parameter f_ij varies for different pairs of beads, random and dissipative DPD pairwise forces are the same for all bead pairs.⁶²

In our model, the interaction force parameters between identical pairs of hydrophobic beads O and T were taken to be the same f_OO = f_TT. To induce better demixing of oil molecules from lipid tails and avoid dispersed oils in the ER-bilayer, we used a slightly larger interaction force parameter f_OT > f_OO = f_TT between oil and tail beads than those between identical beads. Throughout the manuscript, we used fixed pairwise interaction force parameters f_ij between bead types i,j as listed in ESI,† Table S1. These molecular interactions are known to represent membrane properties.⁶¹ Our energy scale is the thermal energy k_BT, where k_B is the Boltzmann constant and T = 298 K is the room temperature.

The basic time scale τ of the model was taken to be 1 ns based on the experimentally measured properties of lipid bilayers, see ref. 60 and the references therein. All DPD simulations were performed using a constant time step of Δt = 0.01τ = 0.01 ns.

In addition to non-bonded interactions between the beads, we introduced intermolecular bonding and bending forces within the lipid and oil chains. A harmonic bond potential

was applied between the adjacent beads of oil and lipid chains at a distance r with an equilibrium bond length r_eq = 0.5d and a bond coefficient k_r = 128k_BT/d². Moreover, a bending potential

V_bend(θ) = k_θ[1 − cos(θ − θ_eq)] (3)

with bending constant k_θ = 15k_BT and equilibrium tilt angle θ_eq = 0, acts on all there consecutive beads in each hydrophobic lipid tail and oil chain.

For all simulations performed here, DPD beads were placed in periodic simulation boxes whose dimensions were chosen in a way to give a bead density 3/d³ = 5.86/nm³. This density is shown to closely estimate water compressibility at room temperature when the interaction parameter f_WW = 25k_BT/d is applied between water beads. We performed three main classes of simulations: spontaneous formation of LDs from initial oil slabs with various values of area-per-lipid, morphological transition of LD lenses to spherical LDs, and LD budding toward the cytosol by inter-monolayer lipid exchange.

Results and discussion

LD structural biogenesis usually occurs in four successive steps: LD nucleation or coalescence i.e. forming nascent symmetrical LD lenses from initially dispersed neutral lipids in the ER-bilayer; LD shape transition from nascent LD lenses to symmetrical spherical LDs; LD budding into the cytosol; and LD fission or detachment from the ER. Our focus, here, is on the two intermediate steps of LD shape transition and budding, before which we identify relevant LD parameters by studying spontaneous LD formation from initial slabs.

LDs, considered here, are either embedded in or attached to the ER membrane and consist of two distinct luminal and cytosolic LD-monolayers, each continuous with its corresponding ER leaflet. Although our membranes are much simpler than actual ER and LD membranes, the lateral bilayer and LD-surrounding monolayers as referred to as the ER-bilayer and LD-monolayers, hereafter. We employed coarse-grained simulations of an LD model to identify LD parameters via spontaneous formation of different LD shapes, and to study LD transition from LD lenses to spherical LDs, and LD budding to the cytosol.

Relevant LD parameters

LD formation, triggered by coalescence of initially dispersed neutral lipids (referred to as oil hereafter) between the ER leaflets, can be understood using the theoretical framework of ternary phase diagram for surfactant-oil–water mixtures.^15,29 Within this framework, LD formation is driven by the reducing contact area between oil and aqueous phases, which lowers the energetic contribution of membrane tensions.

ER-attached LDs consist of two membrane interfaces: the ER-bilayer that separates the two aqueous media—the cytoplasm and the ER lumen; and the two LD-monolayers each separating oil from an aqueous medium. Due to oil hydrophobicity, LD-tensions (at the oil–water interface) are generally larger than the ER-tension (at the water–water interface). As a result, LDs naturally tend to reduce the surface area of their LD-monolayers and expand that of the ER-bilayer. Therefore, the initial state of fully-spread neutral lipids between the ER-leaflets with maximum monolayer area—best represented by a slab structure—spontaneously transforms to an LD lens with smaller monolayers and a larger bilayer. The slab-to-lens transition provides a proper platform to identify relevant LD parameters.

The slab consists of two flat monolayers which are exposed to water from above and below and cover a central oil slab with their hydrophobic chains facing the oil, see Fig. 1(a). Slab monolayers are made from N_lip lipids, equally partitioned between the two monolayers of a fixed width L. The slab-to-lens transition and the resulting LD shapes depend on four parameters: The size of the droplet i.e. the oil volume; the size of the membrane i.e. the area A_s of the slab; and the bilayer and monolayer tensions.

Fig. 1 LD parameters. (a) Initial slab morphology of width L with surrounding water beads (blue). (b) Two slabs of size D = 24.2 nm with different projected areas-per-lipid, A_lp, and corresponding monolayer tensions. (c) LD shapes including slabs, spherical LDs, and LD lenses, spontaneously formed from initial slabs of various sizes with different values of LD parameters, D and A_lp.

The size of an LD is given by the volume V of its oil content, that is the volume of the neutral lipids. In our model, the size of every LD structure, including different LD lenses and the initial slab configuration, with a given oil volume V is defined as the external diameter of a spherical monolayer enclosing the same oil volume, where l_m is the thickness of the LD-monolayer.

The area A_s of the slab is defined as the surface area of one of the equal flat monolayers of the initial slab. The slab area A_s = L², together with the LD size D, determine the surface areas of the ER-bilayer and LD-monolayers of the resulting LDs. For a given D, the minimum slab area, A_sm = πD²/4 + L², is required to form spherical LDs in a box of width L, below which initial slabs retain their flat shape, irrespective of the ER and LD-tensions.

For a given lipid composition, membrane tension has been shown to be directly proportional to its area-per-lipid, for both lipid bilayers and monolayers.^61,63 The area-per-lipid is defined as the average surface area occupied by one lipid molecule. Each slab with a given number of lipids N_lip has a fixed projected area-per-lipid A_lp = 2A_s/N_lip, defined as the average area-per-lipid projected into the horizontal plane. The slab monolayer tension Σ_m is proportional to A_lp, see Fig. 1(b) for small slabs with D = 24.2 nm and A_s = 655.4 nm². Similar behavior was observed for flat bilayers, see ESI,† Fig. S2 for more details. The bilayer and monolayer tensions can be thus replaced by their areas-per-lipid. For a given area A_s, each slab is characterized by two independent parameters: The LD size D, and the projected area-per-lipid A_lp. These parameters determine the resulting LD structure that spontaneously forms from the initial slab.

To see how these parameters affect LD shapes, we simulated slabs of a fixed area A_s = 2304 nm² and different LD sizes D = 29.5, 31, and 33.2 nm. For a given LD size D, the initial flat monolayers with a small N_lip corresponding to a large A_lp, maintained their flat shapes stabilized by large monolayer tensions, see the right snapshots in Fig. 1(c). Monolayers with more lipids i.e. a smaller A_lp, spontaneously transformed into LD lenses as seen in the middle columns in Fig. 1(c). Further increase in the number of lipids of the initial slab eventually led to the spontaneous formation of spherical LDs with smallest A_lp's. LD lenses and spherical LDs that enclose a central axisymmetric oil compartment have been also observed experimentally.^1,16

Similar to previous studies,^11,14,18 our results indicate that the LD shape is related to its size and ER and LD values of area-per-lipid, and thus ER and LD-tensions. Theoretical models have linked the shape of the symmetrical LD lenses to the LD size.¹¹ For a given tension imbalance between the ER and LD membranes, spherical LDs only formed for a relatively large LD volume making LD formation a size-dependent process. However, we observed nano-sized spherical LDs for sufficiently small values of A_lp. Even smaller LDs were detected in multiple spherical structures we obtained for small values of A_lp. Some of these LDs are shown in Fig. S3 (ESI†). We also observed all LD shapes for three different nano-sized droplets we simulated here. Therefore, unlike the conventional view, our results suggest size-independent transition pathways of different character between the LD shapes.

To better understand LD formation, we explored the hypothesis that symmetrical spherical LDs can form via a size-independent continuous transition from LD lenses driven by energetic forces of ER and LD-tensions. Starting from an LD lens of size D = 29.5 nm, spontaneously formed from a slab, we changed the box width L to study LD shape transition in a systematic manner. By calculating membrane tensions and the corresponding A_l's in the ER-bilayer and LD-monolayers, we explored the transition from nascent LD lenses to spherical LDs and revealed how the subtle interplay between the ER and LD areas-per-lipid govern this transition and lead to tense LD-monolayers and nearly tensionless ER-bilayer. See Methods and ESI† for more details on the simulation protocols.

Tension-induced transition of LD lenses to spherical LDs

Each LD lens is composed of two symmetrical monolayers connected to the lateral ER-bilayer. The shape of an LD lens of size D is characterized by the lens angle θ obtained from fitting two spherical caps to its cytosolic and luminal monolayers, see Fig. 2(a). We obtained different LD lenses of a fixed size D = 29.5 nm with lens angles in the range 46° < θ < 86°, each stabilized for a given box width L with the corresponding A_lp, implying the stability of intermediate lenses, see Fig. 2(b).

Fig. 2 Tension-induced transition of LD lenses to spherical LDs. (a) A model LD and a schematic symmetrical LD with LD angle θ. (b) Morphological transition of an initial LD lens of size D = 29.5 nm to an almost spherical LD. (c) ER and LD-tensions for the LD of size D = 29.5 nm. (d) Equal tensions in the cytosolic and luminal LD-monolayers verify our accurate tension calculation. Membrane tensions are proportional to A_l for both the ER-bilayer and LD-monolayers.

For each intermediate lens with angle θ, we calculated individual ER and LD-tensions, Σ_b and Σ_m, respectively. The results are shown in Fig. 2(c). The LD-tension was taken to be the average Σ_m = 0.5(Σ_cm + Σ_lm) of its almost equal components in the cytosolic and luminal monolayers of the symmetrical LD, see the top panel in Fig. 2(d). For the wide LD lens with the smallest angle θ = 46°, we found large ER and LD-tensions of Σ_b = 15.5 mN m⁻¹ and Σ_m = 11.7 mN m⁻¹, with the ER-bilayer being more tense than the LD-monolayers. LD lenses in boxes of smaller width L and thus smaller A_lp's exhibited lower ER and LD-tensions. Upon approaching more spherical lenses with larger lens angles, the initially tense ER-bilayer exhibited greater relaxation of tension compared to the curved LD-monolayers resulting in a tension crossover, at about θ = 50° for an LD size D = 29.5 nm, beyond which the ER-tension fell below the LD-tension. The same trend persisted for larger lens angles leading to an almost tensionless ER-bilayer (Σ_b ≈ 0) for nearly spherical LDs with θ > 85°, wherein LD-monolayers maintained a relatively large tension of about 4 mN m⁻¹. To understand this behaviour and reveal its underlying molecular mechanism, we focused on the molecular structures of ER and LD membranes by calculating their actual areas-per-lipid A_l.

The spontaneous formation of LDs from initial slabs was observed to depend on the projected area-per-lipid A_lp, which is related to the ER and LD-tensions of the final LDs. Membrane tensions are, however, proportional to A_l values in the ER-bilayer and the LD-monolayers. Upon decreasing L, the decrease in A_lp is nearly equally shared between the flat ER-bilayer and the curved LD-monolayers as reflected in A_l's of ER and LD membranes, see the bottom panel of Fig. 2(d). Both ER and LD-tensions linearly decrease with reducing areas-per-lipid. The slope of the fitted line is proportional to the bending rigidity of the membrane and is thus nearly two folds larger for the ER-bilayer compared to the LD-monolayers. It is indeed this difference in the slope of the linear relations between A_l's and membrane tensions in the ER-bilayer and LD-monolayers that gives rise to the stronger drop of the ER-tension compared to LD-tension and results in tensionless ER-bilayer and tense LD-monolayers of spherical LDs. The bottom panel of Fig. 2(d) displays the fitted lines to the ER and LD-tensions versus their corresponding A_l for an LD lens of size D = 29.5 nm. As the shape of the LD lens approaches a sphere, A_l decreases about 0.06 nm² in both the ER-bilayer and LD-monolayers. As a result, membrane tension drops more strongly in the ER-bilayer than in the LD-monolayers eventually leading to a nearly tensionless ER with Σ_b = 0.3 mN m⁻¹ and a tense LD with Σ_m = 4.1 mN m⁻¹ for the almost spherical LD with θ = 86°.

In the absence of the ER-bilayer and LD-monolayers at the liquid interfaces, the ternary phase model relates interfacial tensions at oil–water and water–water interfaces from the force balance equation 2Σ_OW cos θ − Σ_WW = 0. In their presence, however, this relation is slightly shifted resulting in a tension difference ΔΣ = 2Σ_m cos θ − Σ_b > 0 for LDs as shown by the black dash-dotted line in Fig. 2(c). This slight deviation is presumably due to the negligible contributions from LD-monolayers bending and the line tension. The tension difference ΔΣ is proportional to the driving force which pulls the ER bilayer and transforms LD lenses to the spherical LDs. Negligible values of ΔΣ indicate that LD transition is primarily regulated by the force balance or mechanical equilibrium which is independent from the LD size. Assuming ΔΣ = 0, the crossover point of the ER and LD-tensions Σ_b = Σ_m is obtained as θ = 60° for all LD sizes. For initial wide lenses with sufficiently small θ, the ER-tension is in the order of 2Σ_m and thus larger than the LD-tension. During the transition from LD lenses to spherical LDs, the luminal and cytosolic LD-tensions remained almost identical, as expected due to the vertical symmetry of the LD lens, indicating the accuracy of our tension calculations, see the top panel in Fig. 2(d).

We have demonstrated quantitatively that the subtle interplay between A_l values in the ER and LD membranes underlies the continuous LD transition from wide LD lenses to spherical LDs via intermediate LD lenses of different widths. This interplay between the ER and LD areas-per-lipid adjusts ER and LD-tensions and fulfils the force balance independent from the LD size. The slight deviation from this balance provides the driving force for the LD transition. While initial wide LD lenses with sufficiently small angles have larger ER than LD-tensions, a subtle sharing of A_l leads to larger LD than ER tension beyond a crossover point, resulting in large LD and negligible ER-tensions for nearly spherical LDs. This trend, observed here for nano-sized LDs, thus appears to be valid for different LD sizes resulting in a continuous, size-independent LD transition. The LD transition thus resembles a second-order phase transition implying the stability of all LD lenses with different values of θ, consistent with experimental observations.^12,14 Our results are consistent with recent theoretical models that explain the transition from LD lenses to spherical LDs in terms of the relative ratio of ER and LD-tensions, where substantially larger ER-tensions—compared to LD-tensions—stabilized LD lenses.^11,15 Contrary to the conventional view, we find that the LD transition is continuous and does not depend on the LD size.

LD shapes from membrane elastic model

Shape transformations of LDs can also be understood on the basis of membrane elasticity theory. Unlike the simplified spherical cap model that assumes LD-monolayers as two spherical caps, see the bottom panel in Fig. 2(a), the elastic model considers them as smooth surfaces whose shapes are found from minimizing the free energy of the ER-LD system. See the ESI† for more details about the elastic model. Here, we combined this model with the spherical cap model to verify our simulation results, in particular the line tensions we computed for LDs. We focused on the same LD of size D = 29.5 nm discussed before. For each LD lens in Fig. 2(b), we used ER and LD-tensions, Σ_b and Σ_m, obtained from molecular simulations, to calculate the line tension λ from the spherical cap model. We then used these tensions to find the corresponding LD shapes from the elastic membrane model. Remarkably, we found a strong agreement between the LD profile (yellow curves) and the outer surface of neutral lipids from the coarse-grained simulations for different LD angles, see Fig. 3. Using the theoretical framework of the elastic membrane model, we thus validated the membrane and line tensions, Σ_b, Σ_m, and λ, obtained from molecular simulations and the spherical cap model. The line tensions we calculated here are also close to those recently found for aqueous droplets in contact with membranes.^57,59

Fig. 3 Comparing LD shapes from coarse-grained simulations and membrane elastic model. Symmetrical LD-lenses of size D = 29.5 nm with different lens angles from molecular simulations, shown in Fig. 2(b), with membrane profiles from elastic membrane model (yellow curves). Membrane profiles are parametrized with the tangent angle ψ = ψ(s) as a function of the arc length s as shown in the top left panel.

LD budding to the cytosol

We have understood, thus far, how initially-nucleated LD lenses continuously transform to spherical LDs. Next, we will explore the possibility of LD budding and emergence to the cytosol by exchanging lipids from the luminal to the cytosolic LD-monolayer and examine how membrane tensions and areas-per-lipid vary during the LD budding. The asymmetry in the lipid number is quantified by the rescaled difference Δ = (N_cm − N_lm)/N_lip between the number of lipids in the cytosolic and luminal LD-monolayers, N_cm and N_lm, respectively. This asymmetry is created by taking lipids from the luminal monolayers into the cytosolic one at constant total number of lipids N_lip. We note that the lipid exchange in nature proceeds via a non-equilibrium process. To calculate the free energies of LD budding, however, we equilibrated each LD with a given Δ.

Starting from the spherical LD of size D = 29.5 nm, we gradually relocated lipid molecules from the luminal monolayer to the cytosolic one, resulting in continuous increase of Δ and subsequent LD budding into the cytosol. Upon increasing Δ, we observed an asymmetric LD structure characterized by two spherical caps as the cytosolic and luminal monolayers with monolayer angles, θ_c and θ_l, respectively. As Δ was increased, θ_l kept decreasing while the cytosolic monolayer bulged out by increasing θ_c until the luminal monolayer became relatively small (Fig. 4(a)), corresponding to the almost fully-budded LD with θ_c = 163°.

Fig. 4 Directional LD budding into the cytosol by inter-monolayer lipid exchange. (a) Schematic LDs budding to the cytosol by lipid exchange with the corresponding cytosolic angle θ_c. (b) ER and LD-tensions as functions of the cytosolic angle θ_c. LD-monolayers have almost equal tensions and ER-tension is negligible during budding. (c) Almost identical A_l's in the two LD-monolayers as a function of the cytosolic angle θ_c. The cytosolic angle θ_c varies almost linearly with the rescaled lipid number asymmetry Δ.

Similar to the LD transition simulations, we computed the ER and LD-tensions for the cytosolic LD budding. Unlike previous experiments that rely on tension imbalance between the two LD-monolayers to explain LD directional budding,¹⁵ we calculated equal tensions of about 4 mN m⁻¹ in the two LD-monolayers, see Fig. 4(b). Interestingly, the ER-bilayer remained almost tensionless during the LD budding from an ER-embedded spherical LD with θ_c = 92° to an almost fully-emerged LD with θ_c = 163°.

Finding almost identical membrane tensions in the cytosolic and luminal LD-monolayers also implies approximately equal A_l values in the two LD-monolayers throughout LD budding. We verified this by computing A_l values for the cytosolic and luminal LD-monolayers as seen in the top panel of Fig. 4(c). We obtained constant and almost identical A_l values in the two LD-monolayers of the emerging spherical LDs. The area of the cytosolic monolayer grows when adding more lipids. The area-per-lipid, however, remains almost constant, presumably with small fluctuations required for budding. Therefore, our results indicate that LD budding does not essentially require a significant tension imbalance between the cytosolic and luminal LD-monolayers. The bottom panel of Fig. 4(c) displays θ_c as a function of Δ for the emerging LDs in Fig. 4(a). An almost linear growth is observed in θ_c upon increasing Δ.

Free energy landscape of LD transition

The underlying mechanism of LD transition can be understood by identifying the free-energy contributions from bending the LD-monolayers, from ER and LD-tensions, and from the line tension. We calculated the free energy using a mean field approach, typically applicable to membranes where large bending and tension energies are minimally affected by thermal fluctuations as observed recently for membrane interactions with aqueous droplets.^57,59 Note that the ER-bilayer, as a flat membrane with no spontaneous curvature, does not contribute to the bending energy. Fig. 5(a) shows the total free energy and its different components for different symmetrical LDs, shown in Fig. 2(b), as a function of θ.

Fig. 5 Free energy landscape of LD transition and budding. (a) Total free energy of an LD of size D = 29.5 nm as the sum of its components as a function of θ during LD transition and θ_c during LD budding. LD budding is energetically neutral as seen from the term μΔ. (b) LD budding model used to describe the budding process of a spherical LD to the cytosol.

During the transition from an LD lens to an almost spherical LD with θ= 86°, the free energy rapidly dropped from about 5 × 10⁴ pN nm to about 1 × 10⁴ pN nm mainly due to the decreasing ER and LD-tensions, which compensated for the small increase in the line energy. In this regard, the LD transition from LD lenses to spherical LDs is similar to the formation of spherical oil droplets in water in the absence of lipid membranes. The latter one is also driven by the surface tension at the oil–water interface. This is reflected in the bending energy curve in Fig. 5(a) which barely contributes to the total free energy of the LD and thus to the LD transition. Overall the free energy of the LD transition indicates the negligible contribution of membrane bending and line tension even for the nano-sized LDs considered here. This is particularly important as the line tensions are generally expected to contribute more significantly at the sub-micron scale.

Free energy landscape of LD budding

LD budding by creating an LD lipid number asymmetry, via lipid exchange at constant N_lip, can be described by a simple budding model. During budding, the area of the cytosolic LD-monolayer increases while that of the luminal monolayer decreases upon exchanging lipids between LD-monolayers. Our simulations demonstrated that this process proceeds in such a way that the A_l values in LD-monolayers and thereby LD-tensions do not change during budding. This suggests that the initial symmetrical spherical LD can bud by retaining its spherical shape as a constant sphere of size R = 0.5D, see Fig. 5(b). As a result, the cytosolic and luminal LD-monolayers can be represented by two spherical caps of this asymmetrical sphere (with respect to the ER-bilayer), each with their monolayer angles, θ_c and θ_l = π − θ_c, respectively. The initial symmetrical sphere is characterized by θ_l = θ_c = π/2. The areas of the cytosolic and luminal LD-monolayers are given by the areas of the spherical caps, A_jm = 2πR²(1 − cos θ_j)(j = c,l).

We learned from free energy calculations of LD transition that energetic contributions of the line tensions are negligible. Therefore, we simplified our calculations by ignoring the line tension for LD budding. The constant bending energy, 8πκ_m, of a spherical LD throughout budding does not contribute to the variation of the free energy. Here, κ_m is the bending stiffness of the LD-monolayers. Moreover, energetic contribution of membrane bending was also demonstrated to be negligible during LD transition. Hence, the free energy of LD budding is composed of the stretching energies of the ER and LD membranes as well as the free energy associated with the lipid exchange from the luminal to the cytosolic LD-monolayer. For the spherical LD of radius R connected to an ER-bilayer of width L, the free energy is given by:

F_b = A_cΣ_b + A_cmΣ_cm + A_lmΣ_lm + μΔ = (L² − πR² sin²θ_c)Σ_b + 2πR²(1 − cos θ_c)Σ_cm + 2πR²(1 − cos θ_l)Σ_lm + μΔ (4)

where Δ is given by Δ = aθ_c + b (fitted to the data in the bottom panel of Fig. 4(c)) and μ is the free energy conjugate of Δ. By setting θ_l = π − θ_c and Σ_lm = Σ_cm, as obtained from simulations (Fig. 4(b)), and minimizing the resulting free energy with respect to θ_c we arrive at:

However, the ER-tension, Σ_b, was demonstrated to be negligible for spherical LDs formed by LD transition (Fig. 2(c)) and remain negligible throughout LD budding (black dots in Fig. 4(b)). The eqn (5) thus implies that μ, and thereby the free energy term μΔ, are negligible. Consequently, LD budding by lipid exchange from the luminal to the cytosolic LD-monolayer, at constant total number of lipids, is energetically neutral.

Fig. 5(a) shows the free energy contributions of membrane tensions, membrane bending, as well as the term μΔ during LD budding for θ_c < 90°. The free energy μΔ, associated to the lipid exchange is plotted with the actual values of μ calculated from eqn (5), see the orange line in Fig. 5(a). The free energy μΔ is negligible compared to the other energetic contributions.

Membrane tensions

The LD-tensions of about 4 mN m⁻¹ that we found here for LDs of about 30 nm in diameter are very close to the LD tensions of ≈2 and ≈4 mN m⁻¹ experimentally measured for LDs from mammalian and Drosophila cells.¹⁴ Comparing two different LD sizes, D = 27.7 and 29.5 nm, we observed almost identical ER and LD-tensions for both LDs, see ESI,† Fig. S5. To examine the size-dependent behavior of LD-tensions, we computed these tensions for a very small LD of size D = 20.3 nm, and a larger one with D = 33.2 nm spontaneously formed from initial slabs with the same method as described earlier. Fig. 6(a) shows ER and LD-tensions of these almost spherical LDs as a function of D. The LD-tensions were observed to slightly reduce with D for small LDs considered here. A more detailed analysis of size-dependent behavior of LD-tensions requires simulations of larger LDs, which will be accessible with future high performance computing.

Fig. 6 Size-dependent ER and LD tensions. (a) Size-dependent membrane tensions in four nearly spherical LDs of different sizes D = 20.3, 27.7, 29.5, and 33.2 nm. LD tension decreases with the LD size D. (b) Line tension along the LD contact line at the interface between neutral lipids, cytosolic water, and luminal water for two LDs of sizes D = 27.7 and 29.5 nm. Negative line tensions for small LD angles approach zero for spherical LDs.

The smallest LD of size D = 20.3 nm, we observed here, confirms our hypothesis on LD formation by size-independent transition process. These small spherical LDs have not been expected to form due to the dominant role of membrane bending¹¹ which is shown to be negligible in our simulations even in the limit of such a small LD size. According to our results, it is the size-independent tension balance that governs formation of different LD shapes through a continuous transition, as predicted earlier by the formation of multi-spherical LDs in Fig. S3 (ESI†).

We also computed the tension along the contact line of the interior oil with cytosolic and luminal aqueous solutions, see the Supplementary text for more details (ESI†). Fig. 6(b) displays the line tensions computed for two LDs of sizes D = 27.7 and 29.5 nm. Large negative line tensions appear to increase with lens angle θ amounting to a final negligible line tension for spherical LDs.

Conclusion

We conducted extensive simulations of a coarse-grained model to gain insight into the mechanisms behind LD formation and budding with a primary focus on the mechanistic aspects of these processes. We began by identifying the relevant LD parameters, via spontaneous formation of LDs from slabs, which paved the way for our hypothesis that symmetrical LDs can form through a tension-induced continuous transition. We simulated several nano-sized LDs of different sizes and obtained various LD shapes for each LD size. The transition between LD lenses and spherical LDs appears to be continuous resembling a second-order phase transition. For small A_lp values together with single spherical LDs, we also identified the presence of smaller multi-spherical LDs, even though they are unstable (Fig. S3, ESI†), suggesting the stability of nano-sized spherical LDs. These findings align with the principles governing oil droplet formation in water in the absence of lipids. This underscores the dominant role of membrane tensions, which are a common factor in both scenarios, and suggests the likely negligible contribution of membrane bending and line tensions in LD formation.

LD formation and budding are known to be regulated by membrane tensions.^11,14,28 Both the formation of symmetrical ER-embedded LDs and the consequent budding of asymmetric LDs depend on ER and LD-tensions. It is shown experimentally that these tensions are proportional to lipid densities of the ER and LD membranes.¹⁸ The interplay between the ER and LD-tensions has been also shown to determine the LD shape and its transition to fully budded LDs that emerge to the cytosol.^11,14 Previous studies have described LD transition and budding in terms of a discontinuous second-order transition which depends on the LD size, membrane tensions, and membrane bending.^11,28

Symmetrical LDs have been suggested to form through a size-dependent process, primarily controlled by membrane bending for small LD sizes and by membrane tensions for larger LDs.¹¹ Consequently, symmetrical spherical LDs only formed for very large LD sizes while small LD sizes resulted in wide LD lenses with small lens angles. In our simulations, however, we observed spherical LDs even for small nano-sized LDs. This observation was further supported by calculating negligible energetic contributions of membrane bending and line tensions, in comparison to the contributions from membrane tensions.

We also demonstrated that the symmetrical LDs can form through a continuous first-order transition at fixed LD size, governed by the mechanical force balance between LD and ER tensions. The slight deviation from this force balance gives rise to a tension difference (ΔΣ > 0), which pulls the ER bilayer and triggers the continuous transition from tense LD lenses to spherical LDs with negligible ER-tensions. Unlike in previous studies where ΔΣ increased with θ,¹¹ our calculations exhibit a reduction in ΔΣ, and thus the pulling force on the ER-bilayer, as we approach the spherical LDs. This makes sense as the initial strong tendency of LD transition is expected to decrease for more spherical LDs, with increasing values of θ.

Similar to LD formation, the subsequent LD budding and directional emergence toward the cytosol has been also linked to the interplay of ER and LD-tensions and the LD size.¹¹ Over a wide range of LD sizes, LD budding is believed to occur through a discontinuous first-order transition that transforms symmetric LD lenses with large ER-tension to asymmetric fully-budded LDs with small ER-tensions. By contrast, we demonstrated here via calculating free energies, that LD budding by lipid exchange between the LD-monolayers is energetically neutral. Unlike LD formation, which is primarily regulated by membrane tensions, LD budding is mainly driven by this lipid exchange, independent from membrane tensions and the LD size. The modest energy required for budding is presumably provided by proteins that act as the driving agents.

Energetically neutral budding by inter-monolayer lipid exchange is primarily due to the identical mechanical tensions in the luminal and cytosolic LD-monolayers. While the number of lipids and the area of the luminal and cytosolic LD-monolayers change by lipid exchange, their lipid densities or areas-per-lipid remain constant. This is the key to constant identical tensions in LD-monolayers and the consequent neutral budding (see eqn (4)) which makes LD budding a tension-independent process.

LD formation through size-independent continuous transitions of LD lenses and energetically neutral LD budding that we have discovered, carry significant implications. Prior understanding suggested that fully-budded asymmetric LDs were the stable structures for sufficiently small ER-tensions across a wide spectrum of LD sizes.¹¹ This understanding made it challenging to explain the stability of experimentally-observed symmetric LDs,^12,14 particularly those with more spherical shapes. However, our new findings, which reveal a continuous size-independent LD transition driven by the reduction in free energy, followed by energetically-neutral budding, provide a clear explanation for the stability of these close-to-spherical LDs.

The role of protein machinery in LD budding remains a contentious topic in the field.¹ While certain studies have asserted that proteins actively participate in LD budding, others have considered proteins as passive observers with no significant role.^1,31,32 The energetically-neutral nature of LD budding that we have uncovered, based on inter-monolayer lipid exchange, provides new insight into the role of proteins in LD budding. On one hand, we demonstrate that LD budding occurs neutrally, with no remarkable energy input required. On the other hand, any energetically-neutral process demands a slight energy input to take place in a quasi-equilibrium manner. Our findings thus suggest that proteins are indeed essential for LD budding, but their energetic contributions are modest. This perspective bridges the gap between the existing opposing views on the role of proteins in LD budding.

Despite the coarse-grained nature of our molecular simulations with soft potentials, we computed very accurate ER and LD-tensions. LD-tensions we computed here for nearly spherical ER-embedded LDs of different sizes, between 3.5 and 5 mN m⁻¹ as seen in Fig. 6(a), are surprisingly close to the LD-tensions of about 3.5 ± 1.5 and 2 ± 0.5 mN m⁻¹ measured experimentally for LDs from mammalian and Drosophila cells, respectively.¹⁴

The negative line tensions, we calculated here for wide LD lenses, increase to approach positive values for almost spherical ER-embedded LDs, see Fig. 6(b). These tensions presumably remain positive during the entire budding of almost spherical LDs. To the best of our knowledge, tensions of contact lines between lipid droplets and lipid bilayers have not been measured experimentally. Whether the negative line tension, reported here for wide LD lenses, is an artefact of the DPD force field or they can actually exist for LDs in the ER awaits further computational and experimental studies. In particular, such future studies could involve other computational modes, such as the Martini force field⁶⁴ among other available coarse-grained models.

Author contributions

A. H. B. designed research; F. K., B. R., and A. H. B. performed research; F. K., A. B., R. G., B. R., and A. H. B. analyzed the data; A. H. B. wrote the paper with substantial contributions from F. K., R. G., and B. R.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

A. H. B. acknowledges support from the Max Planck Society within the framework of Max Planck Partner Group, and from the European Molecular Biology Organization, grant EMBO IG 5032. B. R. acknowledges support from the National Science Centre, Poland, grant no 2020/39/B/NZ1/00377. Open Access funding provided by the Max Planck Society.

Notes and references

1. J. A. Olzmann and P. Carvalho, Nat. Rev. Mol. Cell Biol., 2019, 20, 137–155.
2. A. R. Thiam, R. V. Farese Jr and T. C. Walther, Nat. Rev. Mol. Cell Biol., 2013, 14, 775–786.
3. A. Pol, S. P. Gross and R. G. Parton, J. Cell Biol., 2014, 204, 635–646.
4. S. Martin and R. G. Parton, Nat. Rev. Mol. Cell Biol., 2006, 7, 373–378.
5. T. C. Walther and R. V. Farese Jr, Annu. Rev. Biochem., 2012, 81, 687–714.
6. F. Wilfling, J. T. Haas, T. C. Walther and R. V. Farese Jr, Curr. Opin. Cell Biol., 2014, 29, 39–45.
7. C. Zhang and P. Liu, Protein Cell, 2017, 8, 796–800.
8. N. Krahmer, R. V. Farese Jr and T. C. Walther, EMBO Mol. Med., 2013, 5, 973–983.
9. P. T. Bozza and J. P. Viola, Prostaglandins, Leukotrienes Essent. Fatty Acids, 2010, 82, 243–250.
10. W. M. Henne, M. L. Reese and J. M. Goodman, EMBO J., 2018, 37, e98947.
11. F. Deslandes, A. R. Thiam and L. Forêt, Biophys. J., 2017, 113, 15–18.
12. A. Chorlay, L. Monticelli, J. V. Ferreira, K. B. Mbarek, D. Ajjaji, S. Wang, E. Johnson, R. Beck, M. Omrane and M. Beller, et al. , Dev. Cell, 2019, 50, 25–42.
13. V. Choudhary, N. Ojha, A. Golden and W. A. Prinz, J. Cell Biol., 2015, 211, 261–271.
14. K. B. M’barek, D. Ajjaji, A. Chorlay, S. Vanni, L. Forêt and A. R. Thiam, Dev. Cell, 2017, 41, 591–604.
15. A. Chorlay and A. R. Thiam, Biophys. J., 2018, 114, 631–640.
16. C. L. Jackson, Curr. Opin. Cell Biol., 2019, 59, 88–96.
17. T. C. Walther, J. Chung and R. V. Farese Jr, Annu. Rev. Cell Dev. Biol., 2017, 33, 491–510.
18. A. Chorlay, L. Forêt and A. R. Thiam, Biophys. J., 2021, 120, 5491–5503.
19. N. Jacquier, V. Choudhary, M. Mari, A. Toulmay, F. Reggiori and R. Schneiter, J. Cell Sci., 2011, 124, 2424–2437.
20. S. Cohen, A. M. Valm and J. Lippincott-Schwartz, Curr. Opin. Cell Biol., 2018, 53, 84–91.
21. A. M. Valm, S. Cohen, W. R. Legant, J. Melunis, U. Hershberg, E. Wait, A. R. Cohen, M. W. Davidson, E. Betzig and J. Lippincott-Schwartz, Nature, 2017, 546, 162–167.
22. O. Adeyo, P. J. Horn, S. Lee, D. D. Binns, A. Chandrahas, K. D. Chapman and J. M. Goodman, J. Cell Biol., 2011, 192, 1043–1055.
23. W. Fei, G. Shui, Y. Zhang, N. Krahmer, C. Ferguson, T. S. Kapterian, R. C. Lin, I. W. Dawes, A. J. Brown and P. Li, et al. , PLoS Genet., 2011, 7, e1002201.
24. J. R. Skinner, T. M. Shew, D. M. Schwartz, A. Tzekov, C. M. Lepus, N. A. Abumrad and N. E. Wolins, J. Biol. Chem., 2009, 284, 30941–30948.
25. J. Zanghellini, F. Wodlei and H. Von Grünberg, J. Theor. Biol., 2010, 264, 952–961.
26. V. Choudhary, G. Golani, A. S. Joshi, S. Cottier, R. Schneiter, W. A. Prinz and M. M. Kozlov, Curr. Biol., 2018, 28, 915–926.
27. T. A. Lagace and N. D. Ridgway, Biochim. Biophys. Acta, Mol. Cell Res., 2013, 1833, 2499–2510.
28. M. Wang and X. Yi, Soft Matter, 2021, 17, 5319–5328.
29. A. R. Thiam and L. Forêt, Biochim. Biophys. Acta, Mol. Cell Biol. Lipids, 2016, 1861, 715–722.
30. A. R. Thiam and F. Pincet, PLoS One, 2015, 10, e0133757.
31. R. Yan, H. Qian, I. Lukmantara, M. Gao, X. Du, N. Yan and H. Yang, Dev. Cell, 2018, 47, 248–256.
32. S. Kim, J. Chung, H. Arlt, A. J. Pak, R. V. F. Jnr, T. C. Walther and G. A. Voth, eLife, 2022, 11, e75808.
33. T. Fujimoto and R. G. Parton, Cold Spring Harbor Perspect. Biol., 2011, 3, a004838.
34. G. Van Meer, D. R. Voelker and G. W. Feigenson, Nat. Rev. Mol. Cell Biol., 2008, 9, 112–124.
35. W. R. Bishop and R. M. Bell, Cell, 1985, 42, 51–60.
36. M. Becuwe, L. M. Bond, N. Mejhert, S. Boland, S. D. Elliott, M. Cicconet, X. N. Liu, M. M. Graham, T. C. Walther and R. V. Farese, J. Cell Biol., 2020, 219, e202006111.
37. S. Kim, J. M. Swanson and G. A. Voth, J. Phys. Chem. B, 2022, 126, 2145–2154.
38. V. Zoni, V. Nieto, L. J. Endter, H. J. Risselada, L. Monticelli and S. Vanni, Front. Mol. Biosci., 2019, 6, 124.
39. V. T. Salo, S. Li, H. Vihinen, M. Hölttä-Vuori, A. Szkalisity, P. Horvath, I. Belevich, J. Peränen, C. Thiele and P. Somerharju, et al. , Dev. Cell, 2019, 50, 478–493.
40. W. Pezeshkian, G. Chevrot and H. Khandelia, Chem. Phys. Lipids, 2018, 211, 93–99.
41. H. Khandelia, L. Duelund, K. I. Pakkanen and J. H. Ipsen, PLoS One, 2010, 5, e12811.
42. V. V. Chaban and H. Khandelia, J. Phys. Chem. B, 2014, 118, 11145–11151.
43. A. Bacle, R. Gautier, C. L. Jackson, P. F. Fuchs and S. Vanni, Biophys. J., 2017, 112, 1417–1430.
44. A. Amirfazli and A. Neumann, Adv. Colloid Interface Sci., 2004, 110, 121–141.
45. R. D. Groot and K. Rabone, Biophys. J., 2001, 81, 725–736.
46. M. Venturoli, M. M. Sperotto, M. Kranenburg and B. Smit, Phys. Rep., 2006, 437, 1–54.
47. T. Auth, D. A. Fedosov and G. Gompper, The Giant Vesicle Book, CRC Press, 2019, pp.169–193.
48. V. Ortiz, S. O. Nielsen, D. E. Discher, M. L. Klein, R. Lipowsky and J. Shillcock, J. Phys. Chem. B, 2005, 109, 17708–17714.
49. J. C. Shillcock and R. Lipowsky, Nat. Mater., 2005, 4, 225–228.
50. W. D. Bennett and D. P. Tieleman, Biochim. Biophys. Acta, Biomembr., 2013, 1828, 1765–1776.
51. K. Müller, D. A. Fedosov and G. Gompper, J. Comput. Phys., 2015, 281, 301–315.
52. F. de Meyer and B. Smit, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 3654–3658.
53. M. Kranenburg, J.-P. Nicolas and B. Smit, Phys. Chem. Chem. Phys., 2004, 6, 4142–4151.
54. H. Noguchi, N. Kikuchi and G. Gompper, Europhys. Lett., 2007, 78, 10005.
55. N. Ramakrishnan, R. P. Bradley, R. W. Tourdot and R. Radhakrishnan, J. Phys.: Condens. Matter, 2018, 30, 273001.
56. M. M. Anila, R. Ghosh and B. Różycki, Soft Matter, 2023, 19, 3723–3732.
57. V. Satarifard, A. Grafmueller and R. Lipowsky, ACS Nano, 2018, 12, 12424–12435.
58. R. Ghosh, V. Satarifard, A. Grafmueller and R. Lipowsky, Nano Lett., 2019, 19, 7703–7711.
59. F. K. Sabet, A. Bahrami and A. H. Bahrami, RSC Adv., 2022, 12, 32035–32045.
60. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
61. B. Różycki and R. Lipowsky, J. Chem. Phys., 2015, 142, 054101.
62. R. D. Groot and P. B. Warren, J. Chem. Phys., 1997, 107, 4423–4435.
63. J. C. Shillcock and R. Lipowsky, J. Chem. Phys., 2002, 117, 5048–5061.
64. S. J. Marrink, H. J. Risselada, S. Yefimov, D. P. Tieleman and A. H. De Vries, J. Phys. Chem. B, 2007, 111, 7812–7824.

---

Footnote

† Electronic supplementary information (ESI) available: Nine supplementary texts, five supplementary figures and seven supplementary tables. See DOI: https://doi.org/10.1039/d3sm01438j

---