Large stress asymmetries of lipid bilayers and nanovesicles generate lipid flip-flops and bilayer instabilities

Abstract

Much effort has been devoted to lipid bilayers and nanovesicles with a compositional asymmetry between the two leaflets of the bilayer membranes. Here, we address another fundamental asymmetry related to lipid densities and membrane tensions. To avoid membrane rupture, the osmotic conditions must be adjusted in such a way that the bilayer membranes are subject to a relatively low bilayer tension. However, even for vanishing bilayer tension, the individual leaflets can still experience significant leaflet tensions if one leaflet is stretched whereas the other leaflet is compressed. Such a stress asymmetry between the two leaflets can be directly controlled in molecular dynamics simulations by the initial assembly of the lipid bilayers. This stress asymmetry is varied here over a wide range to determine the stability and instability regimes of the asymmetric bilayers. The stability regime shrinks with decreasing size and increasing membrane curvature of the nanovesicle. In the instability regimes, the lipids undergo stress-induced flip-flops with a flip-flop rate that increases with increasing stress asymmetry. The onset of flip-flops can be characterized by a cumulative distribution function that is well-fitted by an exponential function for planar bilayers but has a sigmoidal shape for nanovesicles. In addition, the bilayer membranes form transient non-bilayer structures that relax back towards ordered bilayers with a reduced stress asymmetry. Our study reveals intrinsic limits for the possible magnitude of the transbilayer stress asymmetry and shows that the leaflet tensions represent key parameters for the flip-flop rates.

---

1 Introduction

Lipid bilayers are built up from two leaflets of lipid molecules. These bilayers have a thickness of only 4 to 5 nanometers but form closed vesicles with much larger sizes, varying from tens of nanometers for nanovesicles to tens of micrometers for giant vesicles. Giant vesicles have become versatile research tools for basic membrane science, bioengineering, and synthetic biology.^1–3 Nanovesicles provide delivery platforms for pharmaceutics, nutrition, as well as cosmetics and are now used as nanocarriers for the delivery of dermal and transdermal drugs,^4,5 of small interfering RNA therapeutics,^6,7 and of mRNA vaccines.^7,8

In addition, lipid bilayers provide the basic building blocks for all cellular membranes including the outer plasma membrane of all cells, the intracellular membranes within eukaryotic cells, as well as extracellular vesicles such as exosomes and microvesicles. Extracellular vesicles are shed off from almost every type of cell⁹ and are crucial for the chemical communication between cells. Exosomes are released from multivesicular bodies within eukaryotic cells, whereas microvesicles bud off from the outer plasma membranes. Extracellular vesicles have been intensely studied as possible biomarkers for diseases^10,11 and as potential drug delivery systems.^12,13

The two leaflets of cellular membranes typically differ in their lipid composition.¹⁴ Much effort has been devoted to mimick this compositional asymmetry in synthetic vesicles.^3,15,16 For giant vesicles, this asymmetry can now be controlled in a systematic and quantitative manner by low densities of membrane-bound proteins.¹⁷ For synthetic nanovesicles, several protocols have been developed to produce such vesicles with two bilayer leaflets that differ in their phospholipid composition.^16,18–23

Here, we address another fundamental transbilayer asymmetry related to lipid density and membrane tension. Each leaflet of a bilayer membrane experiences a certain leaflet tension which determines how much this leaflet is stretched or compressed. In general, the two leaflet tensions are different and their difference defines the stress asymmetry of the bilayer. So far, stress asymmetries have been largely ignored in experimental studies even though they are relevant for all bilayer membranes. On the other hand, stress asymmetries can be systematically studied in molecular dynamics simulations. Indeed, a few consequences of small leaflet tensions and stress asymmetries have been studied in previous simulations.^24–27 Here, we consider an extended range of leaflet tensions and show that large tensions and asymmetries lead to stress-induced flip-flops of lipids from one leaflet to the other and to the formation of transient nonbilayer structures as illustrated in Fig. 1. As described further below, the nonbilayer structures subsequently undergo a self-healing process, which acts to restore the ordered bilayer structure with a lower stress asymmetry. In this way, our study reveals intrinsic stability limits of bilayer membranes which become unstable when the stress asymmetry exceeds certain limiting values.

Fig. 1 Transient nonbilayer structures arising from instabilities of tensionless bilayers and nanovesicles: (a) globular micelle expelled from the highly compressed upper leaflet of a planar bilayer; and (b) cylindrical micelle expelled from the highly compressed outer leaflet of a nanovesicle. These transient conformations decay via lipid flip-flops which lead to a self-healing process that acts to restore the bilayer structure, see later figures further below as well as the time-lapse Movies M1–M3 in the ESI.† The color coding distinguishes the red-green lipids, initially assembled in the upper and outer leaflets, from the purple-blue lipids, initially assembled in the lower and inner leaflets, which facilitates the visualization of flip-flops.

It is important to distinguish the individual leaflet tensions from the bilayer tension which is equal to the sum of the two leaflet tensions. To avoid membrane rupture, the osmotic conditions have to be adjusted in such a way that both cellular and synthetic membranes experience a relatively low bilayer tension. However, even for vanishing bilayer tension, the individual leaflets can still experience significant leaflet tensions if one leaflet is stretched by a positive leaflet tension whereas the other leaflet is compressed by a negative leaflet tension.

In molecular simulations, the stress asymmetry between the two bilayer leaflets can be controlled by the number of lipids that are initially assembled within the leaflets. Both leaflets are tensionless for a unique reference state in which each leaflet is built up from a specific number of lipids. For planar bilayers, this reference state is easy to determine because it is provided by two leaflets with the same number of lipids. For nanovesicles, the reference state is more difficult to identify because it consists of two leaflets with different lipid numbers. Furthermore, the stress-induced flip-flops and bilayer instabilities are shown to depend on membrane size and curvature. To address this dependence, we study both planar bilayers with vanishing mean curvature and spherical nanovesicles with two different diameters and mean curvatures.

When a lipid molecule flips from one leaflet to the other, the hydrophilic head group of the lipid has to move across the hydrophobic core of the bilayer. The corresponding free energy barrier determines the flip-flop rate via an Arrhenius-like law.^28,29 The flip-flops change the stress asymmetry between the two leaflets of a bilayer. Indeed, when we start from a tensionless bilayer with one compressed and one stretched leaflet, flip-flops from the compressed to the stretched leaflet act to reduce the stress asymmetry. Vice versa, the stress asymmetry can be increased by reshuffling lipids from one leaflet to the other.

From an experimental point of view, the measurement of flip-flop rates has been quite challenging.¹⁵ For phospholipids, a wide range of flip-flop rates has been reported, with the typical time scale for one flip-flop varying from minutes to days.^15,30–38 Likewise, the published values for the inverse flip-flop rates of cholesterol and related sterols were observed to vary from seconds^39,40 to milliseconds.⁴¹ One important conclusion from our simulation study is that the stress asymmetry between the two leaflet tensions of a bilayer membrane represent key parameters for the flip-flop rates. In the experimental studies, these leaflet tensions are difficult to control and play the role of ‘hidden variables’ that contribute to the uncertainty about the experimentally measured flip-flop rates.

Our paper is organized as follows. First, we study the behavior of planar bilayer patches with a constant total number of lipids. We focus on tensionless bilayers and vary the leaflet tensions by reshuffling lipids from one leaflet to the other. For these bilayers, we determine the tension range for which the bilayer remain (meta)stable on the time scale of the simulations. Outside of this tension range, the lipids undergo flip-flops and the bilayers exhibit structural instabilities. To characterize the tension-dependence of the flip-flop rate, we determine the cumulative distribution functions for the statistics of the first flip-flop event. For planar bilayers, this distribution is well fitted by an exponential distribution, which depends only on one fit parameter, the flip-flop rate. We then discuss the structural instabilities of the bilayer and the resulting transient conformations as in Fig. 1a. In the remaining sections, we extend our computational approach to spherical nanovesicles with two different diameters, for which we find qualitatively similar behavior as for planar bilayers as illustrated in Fig. 1b. One notable difference is the statistics for the onset of flip-flops. The corresponding cumulative distribution has a sigmoidal shape which implies that the flip-flop rate increases with the age of the metastable bilayer. At the end, we summarize our results and discuss extensions to bilayers with multiple lipid components and to cellular membranes with membrane proteins that act as lipid pumps.

2 Results and discussion

2.1 Stability of planar bilayers

First, let us consider planar bilayers with vanishing mean curvature. To avoid bilayer poration and membrane rupture, each bilayer is kept at very low bilayer tension by adjusting the total number of lipids assembled in the bilayer to the size of the simulation box. We focus on one-component bilayers with a total number of 1682 lipids assembled in both leaflets. Such a bilayer is symmetric when each leaflet consists of 841 lipids but becomes asymmetric when the two leaflets contain different number of lipids. When the upper leaflet contains N_ul = 870 lipids and the lower leaflet N_ll = 812 lipids, for example, the upper leaflet is compressed and experiences a negative leaflet tension Σ_ul < 0 whereas the lower leaflet is stretched and subject to a positive leaflet tension Σ_ll > 0, see upper panel of Fig. 2. The computation of the leaflet tensions is based on the stress profile across the bilayer as described in the Methods section and illustrated in Fig. S1 of the ESI.†

Fig. 2 Stable and unstable states of planar bilayers, assembled from N_ul lipids in the upper leaflet and N_ll = 1682 − N_ul lipids in the lower leaflet. The blue and red data represent the upper and lower leaflet tensions, Σ_ul and Σ_ll, as functions (a) of N_ul and (b) of area per lipid, a_ul. The green data correspond to the bilayer tension Σ = Σ_ul + Σ_ll, which is close to zero. During the run time of 12.5 μs, we observed no flip-flops within the stability regime (white), corresponding to 737 ≤ N_ul ≤ 945. The left vertical line at N_ul = 957 represents the instability line at which the bilayers become unstable via flip-flops from the compressed upper to the stretched lower leaflet. The right vertical line at N_ul = 725 corresponds to the instability line at which the bilayers become unstable via flip-flops from the compressed lower to the stretched upper leaflet. Each data point represents the mean value of three independent simulations. Apart from the left-most data, the error bars are smaller than the size of the symbols used for the data points. The numerical values of the data are given in Table S1 (ESI†).

In Fig. 2a, the leaflet tensions Σ_ul and Σ_ll are plotted as a function of the lipid number N_ul in the upper leaflet. In Fig. 2b, the same data are shown as a function of the area per lipid, a_ul, in the upper leaflet. For both leaflets, the numerical values of lipid numbers, areas per lipid, and leaflet tensions are provided in Table S1 (ESI†). The areas per lipid are in units of d², the tensions in units of k_BT/d² where d is the diameter of the beads used in our coarse-grained simulation approach and k_BT represents the thermal energy, see Methods section. To simplify the mathematical formulas, we will use the short-hand notation γ ≡ k_BT/d² for the basic tension unit.

During the simulation time of 12.5 μs, the lipids did not undergo any flip-flops in the stability regime which is characterized by lipid numbers N_ul within the range 945 ≥ N_ul ≥ 737 and by upper leaflet tensions Σ_ul with −1.86γ ≤ Σ_ul ≤ + 2.00γ. Flip-flops from the upper to the lower leaflet were observed, however, for N_ul ≥ 957, corresponding to the left instability regime in Fig. 2. For N_ul = 957, which defines the left instability line in Fig. 2, the upper leaflet was initially compressed with Σ_ul = −2.21γ, which induced flip-flops from the compressed upper to the stretched lower leaflet. As the lipid number in the upper leaflet was increased to N_ul > 957, that is, into the left instability regime in Fig. 2, flip-flops became more frequent. On the other hand, for N_ul ≤ 725 and Σ_ul ≥ 2.13γ, which defines the right instability region in Fig. 2, lipids underwent flip-flops from the compressed lower leaflet to the stretched upper one.

2.2 Tensionless and asymmetric bilayers

Inspection of Fig. 2 shows that both leaflet tensions vanish for the symmetric bilayer with N_ul = N_ll = 841 ≡ N₀ and the same lipid areas a_ul = a_ll = 1.271 ≡ a₀ in both leaflets. For this tensionless reference state, both leaflets have the same area A₀ = N₀a₀ = 1023.5d². As we reshuffle lipids from one leaflet to the other, the leaflet areas A_ul = N_ula_ul and A_ll = N_lla_ll remain essentially constant and close to A₀ = N₀a₀. This feature reflects the constraint of vanishing bilayer tension, Σ = Σ_ul + Σ_ll = 0, during the reshuffling of the lipid molecules between the two leaflets and implies the relation a_ll = a_ulN_ul/N_ll between the lipid areas in the two leaflets. The tensionless bilayer becomes asymmetric when the two leaflets contain different lipid numbers N_ul and N_ll ≠ N_ul. This lipid number asymmetry is equivalent to a density asymmetry as described by different areas per lipid in the two leaflets. In addition, the two leaflet tensions Σ_ul and Σ_ll = −Σ_ul have the same magnitude but opposite signs, corresponding to one stretched and one compressed leaflet. In what follows, we will focus on the case of a compressed upper leaflet and a stretched lower one, corresponding to N_ul > N₀ and N_ll = 1682 − N_ul < N₀. In this case, the upper leaflet tension Σ_ul is negative and the lower leaflet tension Σ_ll is positive.

The stress-asymmetry of a planar bilayer can be described by the difference in the leaflet tensions as given by

ΔΣ_pl ≡ Σ_ul − Σ_ll. (1)

For a tensionless bilayer with Σ_ll + Σ_ul = 0 as considered here, the stress asymmetry ΔΣ_pl is equal to 2Σ_ul = −2Σ_ll. Another quantitative measure for this asymmetry is provided by the first moment of the stress profile as defined by eqn (10) in the Methods section.

2.3 Flip-flops in planar bilayers

To determine the kinetic rates of lipid flip-flops, we considered several ensembles of bilayers. Each ensemble consisted of more than 120 bilayers that were assembled from the same lipid numbers N_ul and N_ll and, thus, experienced the same leaflet tensions as long as they remained in their initial metastable states. More precisely, the bilayers experienced the same initial leaflet tensions until time t₁, at which the first flip-flop occurred. The statistics of t₁ can be described by the cumulative distribution function P_cdf(t) that represents the probability that the first flip-flop occurs at time t₁ ≤ t. From a numerical point of view, the cumulative distribution function P_cdf(t) is more reliable than the probability density function dP_exp/dt because, in contrast to the latter function, P_cdf(t) does not require any binning of the data. Furthermore, in the present context, the complementary probability distribution 1 − P_cdf(t) represents the survival probability of the metastable bilayer state without any flip-flops up to time t.

The measured distributions P_cdf(t) are displayed in Fig. 3 for planar bilayers with N_ul = 986, 1015, and 1073 lipids, corresponding to the black circles, red squares, and blue diamonds, respectively. All three data sets can be well fitted by an exponential distribution of the form

P_exp(t) = 1 − exp(−ω_plt) (2)

which involves only a single fit parameter, the flip-flop rate ω_pl for planar bilayers. The inverse rate ω_pl⁻¹ is equal to the average first flip-flop time 〈t₁〉, which also represents the average life time of the metastable bilayer states. Note that the exponential distribution in eqn (2) vanishes for t = 0 and approaches the limiting value P_exp(t) ≈ 1 for large t, as required for any cumulative distribution function, and that the associated probability density function dP_exp/dt is normalized to one.

Fig. 3 Cumulative distribution function P_cdfversus time t, for a planar bilayer with N_ul lipids in the upper leaflet and N_ll = 1682 − N_ul lipids in the lower leaflet. Three sets of data for N_ul = 986 (black circles), N_ul = 1015 (red squares), and N_ul = 1073 (blue diamonds). These data sets are well fitted, using least squares, by an exponential distribution (broken lines) as in eqn (2), which involves only a single fit parameter, the flip-flop rate ω_pl. Each data set represents the outcome of more than 120 statistically independent simulations, see Table S2 (ESI†), which also provides the numerical parameter values of the bilayers. Inset: monotonic increase of the flip-flop rate ω_pl with the absolute value |ΔΣ_pl| = |Σ_ul − Σ_ll| of the stress asymmetry between the two leaflets.

In order to undergo a flip-flop, a lipid has to overcome the free energy barrier arising from the hydrophobic core of the bilayer. When we compress the upper leaflet and stretch the lower one, the flip-flop rate increases monotonically with the stress asymmetry as shown in the inset of Fig. 3. In this inset, the observed flip-flop rate ω_pl in the planar bilayer is displayed versus the absolute value |ΔΣ_pl| of the stress asymmetry as defined by eqn (1); the numerical values of stress asymmetry and flip-flop rate are in Table S2 (ESI†). Thus, our simulations imply that the free energy barrier for flip-flops is decreased as we increase the stress asymmetry |ΔΣ_pl| between the two leaflets of a tensionless bilayer.

2.4 Instability and self-healing of planar bilayers

In addition to the flip-flops of lipid molecules, we also observed structural instabilities of the bilayers as we increased the lipid number N_ul in the upper leaflet beyond the left instability line at N_ul = 957 in Fig. 2. One example for such a structural instability is shown in Fig. 4 and in the time-lapse Movie M1 (ESI†). In this example, the bilayer was initially assembled with N_ul = 986 lipids in the upper leaflet and N_ll = 1682 − N_ul = 696 lipids in the lower leaflet. The structural instability of the bilayer proceeded by the expulsion of about 100 lipids from the compressed upper leaflet, which then formed a globular micelle.

Fig. 4 Structural instability and self-healing of a tensionless planar bilayer. At time t = 0, the bilayer is initially assembled from N_ul = 986 red-green lipids in the compressed upper leaflet and from N_ll = 696 purple-blue lipids in the stretched lower leaflet: (a) at t = 200 ns, the metastable bilayer bulges towards the upper leaflet; (b) at t = 1000 ns, a globular micelle has been formed from about 100 red-green lipids that were expelled from the upper leaflet; (c) at t = 1120 ns, red-green lipids move towards the stretched lower leaflet along the contact line between micelle and bilayer; and (d) this lipid exchange leads to a self-healing process of the bilayer that is completed at t = 1700 ns. At this time point, 93 red-green lipids have moved from the upper to the lower leaflet. The restored bilayer remains stable without flip-flops until the end of the simulations at t = 12.5 μs. For more details on the time course of this process, see time-lapse Movie M1 (ESI†).

The time-lapse Movie M1 (ESI†) can be decomposed into four main stages as shown in Fig. 4. After the initial assembly, the bilayer remains in a metastable state without flip-flops for hundreds of nanoseconds and curves towards the upper leaflet, see the simulation snapshot after 200 ns in Fig. 4a, before it becomes unstable. After 1000 ns, the upper leaflet has expelled some lipids that form a globular micelle, see Fig. 4b. The micelle consists of red-green lipids that were initially assembled in the compressed upper leaflet. Along the contact line between micelle and bilayer segment, we observe frequent exchanges of lipids towards the lower leaflet as can be concluded from Fig. 4c. This lipid exchange increases the number of lipids in the lower leaflet and decreases the number of lipids in the upper leaflet. Finally, after 1700 ns, the conventional bilayer structure has been restored, see Fig. 4d, where the lower leaflet now contains both the blue-purple lipids, from which this leaflet was initially assembled, and 93 red-green lipids that moved from the upper to the lower leaflet. After this self-healing process, the bilayer remained again stable without flip-flops until the end of the simulations at time t = 12.5 μs. Note that the expelled lipids in Fig. 4b form a micelle rather than a bilayer enclosing a certain amount of water, in contrast to the much larger blebs released from the plasma membranes of cells.

In order to study how the lipid flip-flops and bilayer instabilities depend on membrane size and curvature, we extend our simulations to spherical nanovesicles with two different diameters. The larger nanovesicle has a diameter of 23.8d or 19 nm, the smaller one a diameter of 16.3d or 13 nm. The computation of these diameters is described in the Methods section. The mean curvatures of the vesicle membranes are equal to the inverse radii of these vesicles. For both vesicle sizes, we find stable and unstable bilayers of the nanovesicles, in close analogy to the planar bilayers.

2.5 Stability of larger nanovesicles with 19 nm diameter

Next, let us consider the larger nanovesicles with a diameter of 19 nm. The bilayers of the nanovesicles are assembled from a fixed total number of 2875 lipids, with N_il lipids in the inner leaflet and N_ol = 2875 − N_il lipids in the outer one. We use the computational approach described in the Methods section to determine the leaflet tensions Σ_ol and Σ_ilviaeqn (6) and the areas per lipid, a_ol and a_il, of the two leaflets viaeqn (9). The computation of the leaflet tensions is again based on the stress profile s(r) across the bilayer which now depends on the radial coordinate r as illustrated in Fig. S2 (ESI†). The radial coordinate measures the distance from the center of the spherical shapes. We focus again on bilayers with bilayer tension Σ = Σ_ol + Σ_il close to zero.

In Fig. 5a, the two leaflet tensions are plotted as a function of the lipid number N_ol in the outer leaflet; in Fig. 5b, the same tensions are plotted versus the lipid area a_ol in the outer leaflet. Inspection of Fig. 5 and the numerical values in Table S3 (ESI†) shows that both leaflet tensions vanish when the outer leaflet contains between 1893 and 1935 lipids. Linear interpolation then leads to tensionless leaflets for N_ol,0 = 1921 lipids in the outer leaflet and N_il,0 = 954 lipids in the inner one. Thus, in contrast to planar bilayers, nanovesicles with two tensionless leaflets are now characterized by very different lipid numbers in the outer and inner leaflet. Indeed, for the larger nanovesicles considered here, the lipid number in the tensionless outer leaflet is more than twice as large as the one in the tensionless inner leaflet. Likewise, for these tensionless leaflets, the two areas per leaflet are also quite different with a_ol,0 = 1.11d² in the outer leaflet and a_il,0 = 1.45d² in the inner leaflet. Because the area per lipid in the tensionless outer leaflet is smaller than the area per lipid in the tensionless inner leaflet, the latter leaflet is more loosely packed than the tensionless outer leaflet. In fact, inspection of Table S3 (ESI†) reveals that the inequality a_ol < a_il applies to the whole stability regime of the larger nanovesicles which implies that the inner leaflet of such a vesicle is always more loosely packed than the outer leaflet.

Fig. 5 Stable and unstable states of larger nanovesicles, assembled from N_ol lipids in the outer leaflet and N_il = 2875 − N_ol lipids in the inner leaflet. The vesicles have a diameter of 23.8d or 19 nm. The blue and red data represent the outer and inner leaflet tensions, Σ_ol and Σ_il, as functions (a) of the lipid number N_ol and (b) of the area per lipid, a_ol, in the outer leaflet. The green data correspond to the bilayer tension Σ = Σ_ol + Σ_il, which is close to zero. During the run time of 12.5 μs, we observed no flip-flops within the stability regime (white), corresponding to 1775 ≤ N_ol ≤ 2095. The left vertical line at N_ol = 2105 represents the instability line at which the lipids start to undergo flip-flops from the compressed outer to the stretched inner leaflet. The right vertical line at N_ol = 1755 represents the instability line at which the lipids start to undergo flip-flops from the compressed inner to the stretched outer leaflet. Each data point represents the mean value of three independent simulations. The error bars are smaller than the symbols used for the data points. The numerical values of the data are provided in Table S3 (ESI†).

When we start from tensionless leaflets and reshuffle some lipids between these leaflets, we obtain nonzero leaflet tensions Σ_ol and Σ_il = −Σ_ol. The two leaflets form bilayers that remain stable for the lipid number range 1775 ≤ N_ol ≤ 2095 and for the leaflet tension range 1.60γ ≥ Σ_ol ≥ −1.94γ in the outer leaflet, see Fig. 5. Flip-flops from the compressed inner to the stretched outer leaflet are observed for N_ol ≤ 1755 and Σ_ol ≥ + 1.78γ, corresponding to the right shaded region in Fig. 5. Flip-flops from the compressed outer to the stretched inner leaflet occur for N_ol ≥ 2105 and Σ_ol ≤ −1.97γ, which defines the left shaded region in Fig. 5. Within these two instability regimes with N_ol < 1755 and N_ol > 2105, we also observe structural instabilities of the lipid bilayers in addition to the flip-flops.

2.6 Stress-asymmetries in nanovesicles

As shown in Fig. 5, the tensionless bilayer of a nanovesicle may experience significant leaflet tensions Σ_ol and Σ_il = −Σ_ol, corresponding to one compressed and one stretched leaflet. As in the case of planar bilayers, the associated stress asymmetry can be described in two ways. First, by the difference of the two leaflet tensions as given by

ΔΣ_ve ≡ Σ_ol − Σ_il. (3)

Second, the stress asymmetry may also be measured via the first moment [scr T, script letter T]_ve of the stress profile s(r), see eqn (11) in the Methods section.

The stress-asymmetries ΔΣ_ve and [scr T, script letter T]_ve become large and negative as we move away from the tensionless leaflets by increasing the lipid number N_ol in the outer leaflet and decreasing the lipid number N_il in the inner leaflet. Furthermore, these stress-asymmetries become large and positive as we move away from the tensionless leaflets by decreasing the lipid number N_ol in the outer leaflet and increasing the lipid number N_il in the inner leaflet. The numerical values of the stress asymmetries for the larger nanovesicles are included in Table S3 (ESI†).

2.7 Lipid flip-flops in larger nanovesicles

We now focus on the left shaded region in Fig. 5, corresponding to compressed outer leaflets and stretched inner ones. Within this instability region, the lipids undergo flip-flops from the outer to the inner leaflets. To determine the kinetic rates, we examined again different ensembles of nanovesicles. Each ensemble consisted of more than 60 nanovesicles that were initially assembled from the same lipid numbers N_ol and N_il and, thus, experienced the same leaflet tensions as long as they remained in their initial states. More precisely, the bilayers experienced the same initial leaflet tensions until time t₁, at which the first flip-flop occurred. As in the case of planar bilayers, the statistics of t₁ is again described by the cumulative distribution function P_cdf(t) that represents the probability that the first flip-flop occurs at time t₁ ≤ t. This distribution function is depicted in Fig. 6, the numerical values of the data are given in Table S4 (ESI†).

Fig. 6 Cumulative distribution function P_cdfversus time t for nanovesicles with a diameter of 19 nm, assembled from N_ol + N_il = 2875 lipids. Three sets of data are displayed with N_ol = 2105 (black circles), N_ol = 2125 (red squares), and N_ol = 2150 (blue diamonds) lipids in the outer leaflet, which belong to the left instability regime in Fig. 5. The three sets of data are well fitted, using least squares, to Weibull distributions (broken lines) as in eqn (4), which depend on two parameters, the shape parameter k and the rate parameter ω_ve. The parameter values as obtained from the Weibull plots in Fig. S3 are included in Table S4 (ESI†). Each data set represents the outcome of at least 70 statistically independent simulations, see Table S4 (ESI†), which also provides the numerical parameter values of the vesicles. Inset: monotonic increase of the rate parameter ω_ve with the absolute value |ΔΣ_ve| of the stress asymmetry as defined by eqn (3).

The cumulative distribution functions in Fig. 6 have a sigmoidal shape with a point of inflection at intermediate times. The sigmoidal shape is qualitatively different from the exponential function in eqn (2) that was used to fit the onset of flip-flops in planar bilayers. Distribution functions with a sigmoidal shape can be obtained by generalizing the exponential distributions in eqn (2) to Weibull distributions as provided by⁴²

P_Wei(t) = 1 − exp[−(ω_vet)^k]. (4)

This distribution depends on two parameters, the rate parameter ω_ve and the dimensionless shape parameter k > 0. For k = 1, the Weibull distribution in eqn (4) becomes identical to the exponential distribution in eqn (2). For k ≠ 1, the empirical Weibull distribution has been applied to a large variety of different processes.^43,44

The average first flip-flop time 〈t₁〉, which also represents the average life time of the metastable bilayer states, is proportional to the inverse rate parameter ω_ve⁻¹. As explained in the Methods section, the shape parameter k and the rate parameter ω_ve can be estimated from so-called Weibull plots of the observed cumulative distribution function.⁴⁵ The Weibull plots for the larger nanovesicles are displayed in Fig. S3 (ESI†), the numerical values for the shape parameter k and the rate parameter ω_ve are given in Table S6 (ESI†). In all cases, the parameter k is found to be greater than one which corroborates the sigmoidal shape of the distributions and implies that the instantaneous flip-flop rate increases with the age of the metastable state.

In mathematical statistics, the instantaneous rate is known as the hazard rate and equal to the ratio of the probability density function dP_cdf/dt to the survival probability 1 − P_cdf(t).^46,47 The only distribution that leads to a constant and age-independent hazard rate is the exponential distribution. Therefore, the sigmoidal shape of the cumulative distribution functions as observed here for the life time of the metastable bilayer states implies ageing of these nanovesicle states.

2.8 Instability and self-healing of larger nanovesicles

When we move beyond the instability lines in Fig. 5, we observe structural instabilities of the bilayers. We focus on the regimes, in which the outer leaflet is significantly compressed by a negative leaflet tension Σ_ol < 0 and the inner leaflet is significantly stretched by a positive leaflet tension Σ_il > 0. During such a structural instability, the transformation of the nanovesicle proceeds via four stages. One example is displayed in Fig. 7. The details of this shape transformation are displayed in time-lapse Movie M2 (ESI†).

Fig. 7 Structural instability and self-healing process of larger nanovesicle with a diameter of 19 nm. At time t = 0, the bilayer is assembled from N_ol = 2105 and N_il = 770 lipids and the vesicle volume is adjusted in such a way that the bilayer tension is close to zero, which leads to a compressed outer leaflet with negative leaflet tension Σ_ol = −1.97k_BT/d²: (a) at t = 780 ns, the compressed outer leaflet exhibits some kinks; (b) at t = 1720 ns, a cylindrical micelle has been formed from about 180 red-green lipids that were expelled from the outer leaflet; (c) at t = 2160 ns, lipids move towards the stretched inner leaflet along the contact line between micelle and bilayer; and (d) at t = 2710 ns, the self-healing process via stress-induced lipid exchange has been completed and 111 red-green lipids have moved to the inner leaflet. The restored bilayer undergoes no further flip-flops until the end of the simulations at t = 12.5 μs. For more details on the time course of this process, see time-lapse Movie M2 (ESI†). Color coding as in Fig. 1.

In Fig. 7, we see typical snapshots of an individual nanovesicle with a diameter of 19 nm. The bilayer of the vesicle is initially assembled from N_ol = 2105 lipids in the outer and N_il = 770 lipids in the inner leaflet. The outer leaflet is compressed by the negative leaflet tension Σ_ol= −1.97γ, the inner leaflet is stretched by the positive leaflet tension Σ_il = 1.94γ. During the first 1300 ns, the bilayer of the vesicle undergoes shape fluctuations that appear to be ‘normal’ until the outer leaflet starts to develop a protrusion by expelling some lipids. This protrusion grows rapidly into a cylindrical micelle, see Movie M2 (ESI†), which reaches its maximal extension after 1720 ns as shown in Fig. 7b. Somewhat later, some lipids move into the inner leaflet, primarily along the contact line between the micelle and the bilayer, see snapshot at 2160 ns in Fig. 7c. This lipid exchange decreases the number of lipids within the compressed outer leaflet and increases the number of lipids within the stretched inner leaflet until 111 red-green lipids have been moved from the outer to the inner leaflet and the ordered bilayer structure has been restored at 2710 ns as shown in Fig. 7d. After this time point, the restored bilayer undergoes no further flip-flops until the end of our simulations at 12.5 μs. Note that the expelled lipids in Fig. 7b form a micelle rather than a bilayer enclosing a certain amount of water, in contrast to the much larger blebs released from the plasma membranes of cells.

2.9 Stability of smaller nanovesicles with 13 nm diameter

To examine membranes with a higher curvature, we now consider smaller nanovesicles with a diameter of 13 nm. The bilayers of these nanovesicles are assembled from a fixed total number of 1317 lipids, with N_il in the inner leaflet and N_ol = 1317 − N_il in the outer one. Using the same approach as for the larger nanovesicles, we compute the stress profiles s(r), as illustrated in Fig. S4 (ESI†), from which we obtain the leaflet tensions Σ_ol and Σ_il with bilayer tension Σ = Σ_ol + Σ_il close to zero, see Fig. 8 and the numerical values in Table S5 (ESI†). Both leaflet tensions vanish for N_ol,0 = 982 and N_il,0 = 335 as obtained by linear interpolation of the data for N_ol = 977 and N_ol = 987. Note that the tensionless outer leaflet now contains almost three times as many lipids as the tensionless inner leaflet. For these tensionless leaflets, the two areas per leaflet are also quite different with a_ol,0 = 1.09d² and a_il,0 = 1.71d², which implies that the tensionless inner leaflet is again more loosely packed compared to the tensionless outer leaflet. In fact, inspection of Table S5 (ESI†) reveals that the inequality a_ol < a_il applies to the whole stability regime of the smaller nanovesicles which implies that the inner leaflet of such a vesicle is always more loosely packed than the outer leaflet.

Fig. 8 Stable and unstable states of spherical nanovesicles, assembled from 1317 lipids with N_ol lipids in the outer leaflet and N_il = 1317 − N_ol lipids in the inner leaflet. The vesicles have a diameter of 16.3d or 13 nm. The blue and red data represent the outer and inner leaflet tensions, Σ_ol and Σ_il, as functions (a) of the lipid number N_ol and (b) of the area per lipid, a_ol, in the outer leaflet. The green data correspond to the bilayer tension Σ = Σ_ol + Σ_il, which is close to zero. During the run time of 12.5 μs, we observed no flip-flops within the stability regime (white), corresponding to 952 ≤ N_ol ≤ 1037. The left vertical line at N_ol = 1041 represents the instability line at which the bilayers start to undergo flip-flops from the compressed outer to the stretched inner leaflet. The right vertical line at N_ol = 947 represents the instability line at which the bilayers start to undergo flip-flops from the compressed inner to the stretched outer leaflet. Each data point represents the mean value of three independent simulations. The error bars are smaller than the symbols used for the data points. The numerical values of the data are provided in Table S5 (ESI†).

When the lipid numbers assembled in the two leaflets differ from those for the tensionless leaflets, we obtain nonzero leaflet tensions Σ_ol and Σ_il = −Σ_ol as displayed in Fig. 8. Inspection of this figure and the numerical values in Table S5 (ESI†) shows that the bilayers of the small nanovesicles are stable when the lipid number N_ol of the outer leaflet lies within the range 952 ≤ N_ol ≤ 1037 and the outer leaflet tension Σ_ol satisfies 0.72γ ≥ Σ_ol ≥ −1.33γ. Note that this Σ_ol-range is now quite asymmetric with respect to Σ_ol = 0, a feature that is directly visible in Fig. 8. The lipid bilayers of the nanovesicles exhibit flip-flops from the compressed inner to the stretched outer leaflet for N_ol ≤ 947 and positive outer leaflet tension Σ_ol ≥ + 0.79γ, corresponding to the right shaded region in Fig. 8. Flip-flops from the compressed outer to the stretched inner leaflet occur for N_ol ≥ 1041 and negative outer leaflet tension Σ_ol ≤ −1.68γ, which defines the left shaded region in Fig. 8. As we move deeper into the two instability regimes with N_ol < 947 and N_ol > 1041, we observe structural instabilities of the bilayers in addition to the flip-flops.

2.10 Lipid flip-flops in smaller nanovesicles

We focus on the left-shaded region in Fig. 8, in which the lipids undergo stress-induced flip-flops from the compressed outer leaflet to the stretched inner one. The corresponding cumulative distribution functions P_cdf(t) are displayed in Fig. 9. These distributions have again a sigmoidal shape and look qualitatively similar to those in Fig. 6, as obtained for the larger nanovesicles. Furthermore, the three sets of data in Fig. 9 corresponding to three different stress asymmetries are again well fitted by Weibull functions as determined by the Weibull plots in Fig. S5 (ESI†). The corresponding numerical parameter values are provided in Table S6 (ESI†).

Fig. 9 Cumulative distribution function P_cdfversus time t for nanovesicles with a diameter of 13 nm, assembled from a total number of N_il + N_ol = 1317 lipids: three sets of data for lipid numbers N_ol = 1041 (black circles), N_ol = 1046 (red squares), and N_ol = 1051 (blue diamonds) in the outer leaflets. All three data sets belong to the left instability regime in Fig. 8. These data are well fitted, using least squares, to Weibull distributions as in eqn (4) (broken lines), which depend on two parameters, the shape parameter k and the rate parameter ω_ve. The parameter values as obtained from the Weibull plots in Fig. S5 (ESI†) are provided in Table S6 (ESI†). Each data set represents the outcome of more than 60 statistically independent simulations, see Table S6 (ESI†), which also provides the numerical parameter values of the vesicles. Inset: Monotonic increase of the rate parameter ω_ve with the absolute value |ΔΣ_ve| of the stress asymmetry defined by eqn (3).

2.11 Instability and self-healing of smaller nanovesicles

As we move into the instability regimes corresponding to the shaded regions in Fig. 8, the bilayers undergo structural instabilities in addition to the lipid flip-flops. One example is displayed in Fig. S6 and Movie M3 (ESI†). The bilayer of the vesicle is initially assembled from N_ol = 1041 lipids in the outer and N_il = 276 lipids in the inner leaflet. The outer leaflet is compressed by the negative leaflet tension Σ_ol = −1.68γ, the inner leaflet is stretched by the positive leaflet tension Σ_il = 1.71γ. During the first 3400 ns, the bilayer of the vesicle undergoes strong shape fluctuations that have a kinky appearance until the outer leaflet starts to develop a protrusion by expelling some lipids. This protrusion grows into a cylindrical micelle, see Movie M3 (ESI†), which reaches its maximal extension after about 4400 ns as shown in Fig. S6b (ESI†). Somewhat later, some lipids move into the inner leaflet, primarily along the contact line between the micelle and the bilayer, see snapshot at 4800 ns in Fig. S6c (ESI†). This exchange of lipids decreases the number of lipids within the compressed outer leaflet and increases the number of lipids within the stretched inner leaflet. This reshuffling process continues until 57 red-green lipids have been moved from the outer to the inner leaflet and the ordered bilayer structure has been restored after 2710 ns, see Fig. S6d (ESI†). The restored bilayer undergoes no further flip-flops until the end of the simulations at 12.5 μs.

2.12 Comparison between different bilayer systems

We now compare the results for the three bilayer systems as provided by the planar bilayer as well as by the larger and smaller nanovesicles. The three bilayer systems also differ in their mean curvature which vanishes for the planar bilayer and has the values (9.5 nm)⁻¹ and (6.5 nm)⁻¹ for the large and small nanovesicles, corresponding to the inverse radii of these vesicles. All three systems exhibit qualitatively similar but quantitatively different behavior.

Reference states with tensionless leaflets. For each system, the most stable state of the bilayer corresponds to the reference state with two tensionless leaflets. The reference state for a planar bilayer is provided by a symmetric bilayer that contains the same number of lipids in both leaflets, see the intersection point of the upper and lower leaflet tensions in Fig. 2a. In addition, the two tensionless leaflets of a planar bilayer have the same area per lipid. Both symmetries are, however, lost for nanovesicles, for which the reference states consist of two tensionless leaflets that differ both in their number of lipids, N_ol,0 and N_il,0 and in their area per lipid, a_ol,0 and a_il,0. For nanovesicles with a diameter of 19 nm, the tensionless outer leaflet contains N_ol,0 = 1921 lipids with a lipid area of a_ol,0 = 1.11d² whereas the tensionless inner leaflet contains N_il,0 = 954 lipids with a lipid area of a_il,0 = 1.45d², see the intersection point of the two leaflet tensions in Fig. 5a. Likewise, for nanovesicles with a diameter of 13 nm, the tensionless outer leaflet contains N_ol,0 = 1921 lipids with a lipid area of a_ol,0 = 1.09d² whereas the tensionless inner leaflet contains N_il,0 = 335 lipids with a lipid area of a_il,0 = 1.71d², see the intersection point of the two leaflet tensions in Fig. 8a. Thus, in contrast to planar bilayers, the most stable bilayer states of nanovesicles with two tensionless leaflets cannot be obtained by symmetry arguments related to the lipid numbers or the areas per lipid.

Stability regimes and instability lines. A comparison of the three bilayer systems also shows that their stability regimes and instability lines correspond to different leaflet tensions and to different first moments of the stress profiles. In Table 1, we characterize the two instability lines and the stability regimes in terms of the leaflet tensions Σ_ul and Σ_ol experienced by the upper leaflets of the planar bilayers and the outer leaflets of the nanovesicles. Inspection of Table 1 shows that the tension range for the stability regime becomes smaller as we increase the mean curvature of the bilayer system. Indeed, this tension range has an extension of 3.86γ for the planar bilayer, 3.54γ for the large nanovesicles and 2.05γ for the small nanovesicles, with the tension unit γ = k_BT/d².

Table 1 Left and right instability lines and stability regime, as displayed in Fig. 2, 5, and 8, in terms of the upper and outer leaflet tensions Σ_ul and Σ_ol of planar bilayers and nanovesicles, respectively. All tensions are given in units of γ = k_BT/d². The mean curvature vanishes for the planar bilayer and increases from (11.9 nm)⁻¹ for the large nanovesicles to (8.13 nm)⁻¹ for the small ones

	Left instability line	Stability regime	Right instability line
Planar bilayers	Σ ul = −2.21γ	−1.86γ ≤ Σul ≤ + 2.00γ	Σ ul = + 2.13γ
Large nanovesicles	Σ ol = −1.97γ	−1.94γ ≤ Σol ≤ + 1.60γ	Σ ol = + 1.78γ
Small nanovesicles	Σ ol = −1.68γ	−1.33γ ≤ Σol ≤ + 0.72γ	Σ ol = + 0.79γ

---

Another quantity by which one can characterize the instability lines and the stability regimes is the first moment of the stress profiles which provides the torque densities, [scr T, script letter T]_pl and [scr T, script letter T]_ve, as computed viaeqn (10) for the planar bilayers and viaeqn (11) for the nanovesicles. In Table S7 (ESI†), we characterize the stability regime together with the two instability lines in terms of these torque densities, confirming that the stability regime shrinks as we increase the mean curvature of the bilayer system.

Finally, we compare the cumulative distribution functions for the onset of flip-flops in the unstable regimes of the nanovesicles, as displayed in Fig. 6 and 9. We rescale the time t by the rate parameter ω_ve for the nanovesicles, thereby introducing the dimensionless time s = ω_vet. As a result, we obtain the master curve in Fig. S7 (ESI†) which shows that all six data sets for the large and small nanovesicles lie on top of each other, further corroborating the accuracy of the deduced flip-flop rates.

3 Summary and outlook

In summary, we studied the stable and unstable states of three bilayer systems as provided by planar bilayers and nanovesicles with two different diameters. In all three systems, the bilayers were tensionless and characterized by (almost) zero bilayer tension whereas the two leaflets experienced variable leaflet tensions of opposite sign. For each of these systems, we identified a stability regime of low leaflet tensions, in which the bilayers do not undergo flip-flops and retain their ordered bilayer structure during the time scale of our simulations (Fig. 2, 5, and 8). Each of these stability regimes are bounded by two instability lines, at which the bilayers start to undergo stress-induced flip-flops. When the leaflet tensions are increased beyond these instability lines, the bilayers undergo both flip-flops and structural instabilities. Comparison of the results for the three bilayer systems (Table 1) shows that the stability regime is reduced when the bilayer membrane is more strongly curved.

The most stable states within the stability regimes are provided by the reference states with two tensionless leaflets. For planar bilayers, these reference states correspond to symmetric bilayers for which the two leaflets contain the same number of lipids and are characterized by the same area per lipid (Fig. 2 and Table S1, ESI†). Both symmetries are lost for nanovesicles, for which the reference states consist of two tensionless leaflets that differ both in their number of lipids and in their area per lipid (Fig. 5 and Table S3 as well as Fig. 8 and Table S5, ESI†) Thus, in contrast to planar bilayers, the most stable bilayer states of nanovesicles with two tensionless leaflets cannot be obtained by symmetry arguments related to the lipid numbers or the areas per lipid.

In all three bilayer systems, the onset of stress-induced flip-flops can be characterized by a cumulative distribution function P_cdf(t), which is complementary to the survival probability of the metastable states of the bilayers as given by 1 − P_cdf(t). For the planar bilayers, the cumulative distribution function is well fitted by an exponential distribution (Fig. 3). For the nanovesicles, the cumulative distribution function has a sigmoidal shape (Fig. 6 and 9), which implies that the instantaneous flip-flop rate is not constant but increases with the age of the metastable bilayer. In the latter case, the cumulative distribution functions are well fitted by a Weibull distribution as described by eqn (4). The flip-flop rates ω_pl and ω_ve deduced from the distribution functions determine the average time scale for the first stress-induced flip-flop that is observed when the bilayer becomes unstable. For the large and small nanovesicles, all six data sets form a single master curve when plotted against the rescaled time s = ω_vet (Fig. S7, ESI†). This master curve corroborates the reliability of the deduced flip-flop rates.

In addition, all three bilayer systems undergo structural instabilities that proceed via four main stages (Fig. 4, 7, and Fig. S6, ESI†). When the upper leaflet of a planar bilayer or the outer leaflet of a nanovesicle are highly compressed, these leaflets expel many lipids that form a globular (Fig. 4b) or cylindrical micelle (Fig. 7b and Fig. S6b, ESI†). Subsequently, a self-healing process sets in that is mediated via flip-flops towards the highly stretched lower and inner leaflets. A large fraction of these flip-flops is observed along the contact lines between the micelles and the remaining bilayers (Fig. 4c, 7c, and Fig. S6c, ESI†). As a result of these flip-flops, the micelles are reintegrated into the bilayers, which then remain in a metastable state until the end of the simulations (Fig. 4d, 7d, and Fig. S6d, ESI†).

In this study, we focused on one-component bilayers for the sake of computational simplicity but our results can be directly generalized to multi-component bilayers as well, provided the bilayer leaflets are fluid and have a laterally uniform composition. The latter requirement does not apply to coexistence regions of the phase diagram, in which the bilayers undergo lipid phase separation. On the other hand, in the one-phase regions of the phase diagram, the leaflet tensions describe the local stress within the leaflets, irrespective of the number of molecular components as long as the leaflets remain in a fluid state and do not phase separate. Thus, for a multi-component bilayer, the different components will have different areas per lipid but the average molecular area in each leaflet will again determine the corresponding leaflet tension and vice versa. Furthermore, for a two-dimensional leaflet, the associated leaflet tension plays the same role as the pressure in a three-dimensional liquid. Therefore, in the absence of hydrodynamic flows within the bilayer membranes, the leaflet tensions must be uniform along the membranes and will again govern the flip-flop rates, which will now differ for different lipid components.

Multi-component bilayers can undergo lipid phase separation as observed for nanovesicles containing the two phospholipids DPPC and POPC.¹⁸ The inner leaflet of these nanovesicles consisted of a uniform POPC-rich phase whereas the outer leaflet underwent lipid phase separation into a DPPC-rich and a POPC-rich phase. The DPPC-rich phase was not fluid but formed a gel phase. In such a situation, the phase-separated leaflet involves a line tension that arises from the domain boundaries between the DPPC-rich and the POPC-rich phase. Such a line tension would make another contribution to the free energy, of the nanovesicle, in addition to the leaflet tensions, but plays no role for bilayer leaflets with a laterally uniform composition as considered here.

Multi-component bilayers could also be studied in order to measure the cumulative distribution function for the onset of flip-flops experimentally. To prepare a population of vesicles with a significant stress asymmetry, one might use the enzymatic conversion of phospholipid headgroups in the outer leaflet as demonstrated by the conversion of the PS to the PE headgroup²¹ and of the PC headgroup to the PS and PE headgroups.²² After such a vesicle population has been prepared, one would look for flip-flops and measure the corresponding flip-flop times, which are likely to exceed the flip-flop times observed in our simulations by several orders of magnitude. Studying a dilute dispersion of liposomes and adding local lipid probes to their membranes, it may also be possible to distinguish flip-flops in the same vesicle from those in different vesicles. In this way, one could determine the cumulative distribution function P_cdf(t) experimentally, in close analogy to our simulation study. If this distribution has a sigmoidal shape as in Fig. 6 and 9, the flip-flop rate is not constant but depends on the age of the nanovesicles.

Our results are also relevant for the behavior of biological membranes which contain membrane proteins that change the flip-flop rates of the lipids. Three families of such membrane proteins have been identified. One family of membrane proteins, so-called scramblases, act to reduce the free energy barrier for the flip-flops, thereby facilitating bi-directional and reversible flip-flops between the two bilayer leaflets.^14,48–51 In addition, two families of membrane proteins that hydrolyze ATP have also been identified: Flippases actively pump lipids from the exoplasmic to the cytoplasmic leaflet of the membranes whereas floppases hydrolyze ATP to pump lipids in the opposite direction, i.e., from the cytoplasmic to the exoplasmic leaflet.^14,52–54

Because the activity of a scramblase is not coupled to ATP hydrolysis, such a membrane protein acts to reduce the stress asymmetry between the two leaflets. In particular, scramblases should further increase the flip-flop rate from a compressed to a stretched leaflet. On the other hand, both flippases and floppases will increase the stress asymmetry. As an example, consider a nanovesicle with reconstituted flippase proteins that move lipids from the inner to the outer leaflet, thereby compressing the outer leaflet and stretching the inner one. When we reach a sufficiently large stress asymmetry between the two leaflets, the nanovesicle is expected to undergo structural instabilities as displayed in Fig. 7 and Fig. S6 (ESI†). Therefore, using nanovesicles with reconstituted flippases, these structural instabilities should become accessible to systematic experimental studies.

4 Methods

4.1 Molecular modelling and interactions

Our coarse-grained simulation approach is provided by the dissipative particle dynamics (DPD) simulation method.^55,56 In DPD, the lipid and water molecules are modelled as soft beads that interact with pairwise additive forces. These forces conserve momentum locally, which ensures a reliable description of hydrodynamics. The aqueous liquid consists of water (W) beads that have a diameter d of about 0.8 nm and thus represent small volumes of the liquid. The lipids consist of three headgroup (H) beads and two hydrocarbaon chains, each of which is built up from six chain (C) beads, corresponding to the SH lipids used in our previous simulation study.²⁵ All beads have the same diameter d, which provides the basic length scale and the cut-off for the DPD forces; for further details of our method, see ref. 25. The DPD method has been applied to a large variety of membrane topics including diffusion through membranes.⁵⁷ We performed all simulations using the MD code LAMMPS (Large Scale Atomic/Molecular Massively Parallel Simulator).⁵⁸

4.2 Simulation parameters and bilayer geometry

We study planar bilayers consisting of 1682 lipids in total. The reference state with two tensionless leaflets is provided by the symmetric bilayer with 841 lipids in each leaflet. Asymmetric bilayers are obtained if the upper and lower leaflet contain different lipid numbers N_ul and N_ll = 1682 − N_ul, while keeping the total number of lipids constant. The cubic simulation box has the base area A = L_xL_y and the height L_z with L_x = L_y varied between 32.0d and 32.2d and the height L_z varied between 46.9d and 46.3d for fixed volume of the box. The lipid areas in the upper and lower leaflet are then obtained via a_ul = A/N_ul and a_ll = A/N_ll which implies the relationship a_ll = a_ulN_ul/N_ll.

The larger nanovesicles are assembled from a total number of 2875 lipids, the smaller ones from a total number of 1317 lipids. Tensionless bilayers with Σ = 0 are obtained by adjusting the volume of the vesicles. The inner and outer radii R_iH and R_oH of the spherical bilayers are estimated from the two maxima of the density profile ρ_H(r) for the hydrophilic H beads, corresponding to the inner and outer headgroup layers. For the larger nanovesicles, this procedure leads to the inner radius R_iH = 9.2d and to the outer radius R_oH = 14.2d. For the smaller nanovesicles, we obtain the inner radius R_iH = 5.3d and the outer radius R_oH = 10.3d.

4.3 Stress profiles, midsurfaces, and leaflet tensions

The stress profile s(z) across a planar bilayer, which depends on the Cartesian coordinate z, is computed as in ref. 59 and 25. The value z = 0 corresponds to the location of the maximum for the density profile of the hydrophobic C beads, which is taken to define the midplane of the bilayer. The leaflet tensions Σ_ul and Σ_ll of the upper and lower leaflet are given byThe sum Σ_ul + Σ_ll of the two leaflet tensions is equal to the bilayer tension Σ, which implies that Σ_ll = −Σ_ul for a tensionless bilayer with Σ = 0.

The stress profile s(r) across a spherical nanovesicls, which depends on the radial coordinate r, is obtained as described previously.²⁶ The midsurface of the spherical bilayer is placed at r = R_mid, corresponding to the r-value, at which the density profile of the hydrophobic C beads has its maximum. As a result, we obtain the radius R_mid = 11.88d for the midsurface of the larger nanovesicles and the radius R_mid = 8.13d for the midsurface of the smaller nanovesicles.

The leaflet tensions Σ_ol and Σ_il of the outer and inner leaflet are

where R_max is an appropriate upper cutoff for the r-integration. The sum Σ_il + Σ_ol of the two leaflet tensions is equal to the bilayer tension Σ, which implies that Σ_il = −Σ_ol for a tensionless bilayer with Σ = 0. The numerical values for the leaflet tensions obtained in this way are displayed in Tables S1–S6 (ESI†).

4.4 Midsurfaces of leaflets and areas per lipid

In order to determine the areas per lipid in the two leaflets, we first compute the radii R_il and R_ol for the midsurfaces of these leaflets. The inner leaflet is bounded by the inner headgroup layer with radius R_iH and the midsurface of the bilayer with radius R_mid. Thus, the midsurface of the inner leaflet has the radiusLikewise, the outer leaflet is bounded by the midsurface of the bilayer and the outer headgroup layer. Thus, the midsurface of the outer leaflet has the radiusUsing these relationships, we obtain R_il = 10.5d and R_ol = 13.0d for the large nanovesicles as well as R_il = 6.7d and R_ol = 9.2d for the small nanovesicles.

Finally, the areas per lipid, a_il and a_ol, in the inner and outer leaflets are computed by dividing the areas of the midsurfaces of the inner and outer leaflets by the number of lipids assembled within these leaflets, which leads to

As a result, we obtain the numerical values of a_il and a_ol in Tables S1–S6 (ESI†).

4.5 First moments of stress profiles

For a planar and tensionless bilayer, the first moment of the stress profile s(z) can be interpreted as a torque density, [scr T, script letter T]_pl, which is directly related to the product of two curvature-elastic parameters, the bending rigidity κ and the spontaneous curvature m. This relationship has the form^24,25,60The numerical values of [scr T, script letter T]_pl are included in Table S1 (ESI†).

For spherical nanovesicles with a tensionless bilayer, the first moment of the stress profile is related to the bending rigidity κ, the spontaneous curvature m, and the radius R_mid of the spherical mid-surface according to²⁶

The computation of the mid-surface radius R_mid is described in the previous paragraph. The stress profile s(r) decays rapidly to zero as we move away from the bilayer, which implies that the r-values that contribute to the integral in eqn (11) are restricted to the vicinity of the bilayer.

4.6 Weibull plots for the onset of flip-flops

The statistics for the onset of flip-flops is encoded in the cumulative distribution function P_cdf(t) which describes the probability that the first flip-flop occurs at time t₁ ≤ t. The complementary distribution 1 − P_cdf(t) represents the survival probability for the metastable state without any flip-flops up to time t. In order to fit the simulation data to the Weibull distribution in eqn (4), we rewrite the latter equation as

ln(−ln[1 − P_Wei(t)]) = k ln(t) + k ln(ω_ve). (12)

Likewise, the simulation data for the cumulative distribution function P_cdf are transformed into ln(−ln[1 − P_cdf(t)]), and these transformed data are plotted versus ln(t), thereby generating so-called Weibull plots. Least squares fitting of these data to Weibull distributions P_cdf(t) = P_Wei(t) as in eqn (12) then leads to estimates for the shape parameter k and for the rate parameter ω_ve of these distributions. The Weibull plots for the large and small nanovesicles are displayed in Fig. S3 and S6 (ESI†), the numerical values for k and ω_ve are included in Tables S4 for the larger and in Table S6 (ESI†) for the smaller nanovesicles.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This research was conducted within the Max Planck School Matter to Life supported by the German Federal Ministry of Education and Research (BMBF) in collaboration with the Max Planck Society and the Max Planck Institute of Colloids and Interfaces. Open Access funding provided by the Max Planck Society.

References

1. Physics of Biological Membranes, ed. P. Bassereau and P. Sens, Springer Nature, 2018.
2. The Giant Vesicle Book, ed. R. Dimova and C. Marques, Taylor & Francis, 2020.
3. R. Lipowsky and R. Dimova, Soft Matter, 2021, 17, 214–221.
4. S. Rai, V. Pandey and G. Rai, Nano Rev. Exp., 2017, 8, 1325708.
5. T. Limongi, F. Susa, M. Marini, M. Allione, B. Torre, R. Pisano and E. D. Fabrizio, J. Nanomater., 2021, 11, 3391–3429.
6. J. A. Kulkarni, D. Witzigmann, S. Chen, P. R. Cullis and R. V. D. Meel, Acc. Chem. Res., 2019,(52), 2435–2444.
7. Y. Suzuki and H. Ishihara, Drug Metab. Pharmacokinet., 2021, 41, 100424.
8. M. Chen, X. Hu and S. Liu, Macromol. Chem. Phys., 2022, 223, 2100443.
9. E. Woith, G. Fuhrmann and M. F. Melzig, Int. J. Mol. Sci., 2019, 20, 5695.
10. L. Barile and G. Vassalli, Pharmacol. Ther., 2017, 174, 63–78.
11. I. Jalaludin, D. M. Lubman and J. Kim, Mass Spec. Rev., 2021, e21749.
12. I. K. Herrmann, M. J. A. Wood and G. Fuhrmann, Nat. Nanotechnol., 2021, 16, 748–759.
13. L. van der Koog, T. B. Gandek and A. Nagelkerke, Adv. Healthcare Mater., 2022, 11, 2100639.
14. Y. Yang, M. Lee and G. D. Fairn, J. Biol. Chem., 2018, 293, 6230–6240.
15. D. Marquardt, B. Geier and G. Pabst, Membranes, 2015, 5, 180–196.
16. J. R. St. Clair, Q. Wang, G. Li and E. London, The Biophysics of Cell Membranes, Springer Nature Singapore, 2017, vol. 19, pp. 1–27.
17. J. Steinkühler, R. L. Knorr, Z. Zhao, T. Bhatia, S. Bartelt, S. Wegner, R. Dimova and R. Lipowsky, Nature Commun., 2020, 11, 905.
18. F. A. Heberle, D. Marquardt, M. Doktorova, B. Geier, R. F. Standaert, P. Heftberger, B. Kollmitzer, J. D. Nickels, R. A. Dick, G. W. Feigenson, J. Katsaras, E. London and G. Pabst, Langmuir, 2016, 32, 5195–5200.
19. H.-Y. Sun, G. Deng, Y.-W. Jiang, Y. Zhou, J. Xu, F.-G. Wu and Z.-W. Yu, Chem. Commun., 2017, 53, 12762–12765.
20. H.-Y. Guo, H.-Y. Sun, G. Deng, J. Xu, F.-G. Wu and Z.-W. Yu, Langmuir, 2020, 36, 12684–12691.
21. C. Drechsler, M. Markones, J.-Y. Choi, N. Frieling, S. Fiedler, D. R. Voelker, R. Schubert and H. Heerklotz, Biophys. J., 2018, 115, 1509–1517.
22. R. Takaoka, H. Kurosaki, H. Nakao, K. Ikeda and M. Nakano, BBA – Biomembr., 2018, 1860, 245–249.
23. K. Kamiya, T. Osaki and S. Takeuchi, Sens. Actuators, B, 2021, 327, 128917.
24. B. Różycki and R. Lipowsky, J. Chem. Phys., 2015, 142, 054101.
25. A. Sreekumari and R. Lipowsky, J. Chem. Phys., 2018, 149, 084901.
26. R. Ghosh, V. Satarifard, A. Grafmüller and R. Lipowsky, Nano Lett., 2019, 19, 7703–7711.
27. M. Miettinen and R. Lipowsky, Nano Lett., 2019, 19, 5011–5016.
28. A. Imparato, J. C. Shillcock and R. Lipowsky, Eur. Phys. J. E: Soft Matter Biol. Phys., 2003, 11, 21–28.
29. G. Parisio, A. Ferrarini and M. M. Sperotto, Int. J. Adv. Eng. Sci. Appl. Math., 2016, 8, 134–146.
30. R. D. Kornberg and H. M. McConnell, Biochemistry, 1971, 10, 1111–1120.
31. B. W. van der Meer, Biomembranes: Physical Aspects, VCH Verlagsgesellschaft, 1993, pp. 97–158.
32. A. V. Victorov, N. Janes, G. Moehren, E. Rubin, T. F. Taraschi and J. B. Hoek, J. Am. Chem. Soc., 1994, 116, 4050–4052.
33. J. Liu and J. C. Conboy, Biophys. J., 2005, 89, 2522–2532.
34. M. Nakano, M. Fukuda, T. Kudo, N. Matsuzaki, T. Azuma, K. Sekine, H. Endo and T. Handa, J. Phys. Chem. B, 2009, 113, 6745–6748.
35. J. Liu, K. L. Brown and J. C. Conboy, Faraday Discuss., 2013, 161, 45–61.
36. D. Marquardt, F. A. Heberle, T. Miti, B. Eicher, E. London, J. Katsaras and G. Pabst, Langmuir, 2017, 33, 3731–3741.
37. M. M. Sperotto and A. Ferrarini, The Biophysics of Cell Membranes, 2017, 29–60.
38. L. Porcar and Y. Gerelli, Soft Matter, 2020, 16, 7696–7703.
39. T. L. Steck, J. Ye and Y. Lange, Biophys. J., 2002, 83, 2118–2125.
40. J. A. Hamilton, Curr. Opin. Lipidol., 2003, 14, 263–271.
41. R. J. Bruckner, S. S. Mansy, A. Ricardo, L. Mahadevan and J. W. Szostak, Biophys. J., 2009, 97, 3113–3122.
42. W. Weibull, J. Appl. Mech., 1951, 18, 293–297.
43. H. Rinne, The Weibull Distribution – A Handbook, Taylor & Francis, 2008.
44. C.-D. Lai, Generalized Weibull Distributions, Springer, 2014.
45. C. Corsaro, G. Neri, A. M. Mezzasalma and E. Fazio, Polymers, 2021, 13, 2897.
46. D. R. Cox, Renewal Theory, Methuen and Co Ltd, 1962.
47. H. M. Taylor and S. Karlin, An Introduction to Stochastic Modelling, Academic Press, San Diego, 1998.
48. P. Comfurius, P. Williamson, E. F. Smeets, R. A. Schlegel, E. M. Bevers and R. F. A. Zwaal, J. Biol. Chem., 1996, 35, 7631–7634.
49. X. Buton, G. Morrot, P. Fellmann and M. Seigneuret, J. Biol. Chem., 1996, 271, 6651–6657.
50. J. Suzuki, M. Umeda, P. J. Sims and S. Nagata, Nature, 2010, 468, 834–840.
51. E. M. Bevers and P. L. Williamson, Physiol. Rev., 2016, 96, 605–645.
52. P. F. Devaux, FEBS Lett., 1988, 234, 8–12.
53. T. T. Sebastian, R. D. Baldridge, P. Xu and T. R. Graham, Biochim. Biophys. Acta, 2012, 1821, 1068–1077.
54. H. M. Hankins, R. D. Baldridge, P. Xu and T. R. Graham, Traffic, 2015, 16, 35–47.
55. R. D. Groot and P. B. Warren, J. Chem. Phys., 1997, 107, 4423–4435.
56. J. C. Shillcock and R. Lipowsky, J. Chem. Phys., 2002, 117, 5048–5061.
57. Z. Xu, L. Gao, P. Chen and L.-T. Yan, Soft Matter, 2020, 16, 3869–3881.
58. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
59. R. Goetz and R. Lipowsky, J. Chem. Phys., 1998, 108, 7397–7409.
60. W. Helfrich, Physics of Defects, North-Holland Publishing Company, Amsterdam, 1981, pp. 715–755.

---

Footnote

† Electronic supplementary information (ESI) available: Supporting Fig. S1–S7, supporting Tables S1–S7, and three Movies M1–M3. See DOI: https://doi.org/10.1039/d2sm00618a

---