Balancing stability and flexibility when reshaping archaeal membranes

To study bolalipid membranes and compare them to bilayers, we extended the Cooke and Deserno model for bilayer membranes (Cooke and Deserno, 2005). In the original model for the bilayer a single bilayer lipid is represented by a chain of three nearly equally sized beads of diameter of $\sim 1\sigma$, where $\sigma$ is our distance unit and roughly maps to $1\mathrm{n}\mathrm{m}$ (Figure 1C left); one bead stands for the head group (cyan), while the others represent the hydrophobic tail (blue). Each adjacent pair of beads in a lipid is linked by a finite extensible nonlinear elastic (FENE) bond. The angle formed by the chain of three beads is kept near 180° via an angular potential with strength $k_0$, instead of the approximation by a bond between end beads of the original model (Cooke and Deserno, 2005).

While lipid heads interact exclusively through volume exclusion, the beads of lipid tails interact via a soft attractive potential of the strength ${ϵ}_{\mathrm{p}}$ and range $\omega$ (Figure 1C left), effectively modelling hydrophobic interaction in an implicit solvent. This interaction strength governs the membrane phase behaviour and can be interpreted as the effective temperature or reduced temperature $T_{\mathrm{e}\mathrm{f}\mathrm{f}}=k_{\mathrm{B}}T/{ϵ}_{\mathrm{p}}$. As the distinction between scaling interactions ($T_{\mathrm{e}\mathrm{f}\mathrm{f}}$) or temperature ($T$) is not important for our analysis (see Appendix 1–Section 14), for simplicity, we refer to $T_{\mathrm{e}\mathrm{f}\mathrm{f}}$ as temperature in the following.

To model a bolalipid molecule, we joined two bilayer lipids so that a lipid molecule is formed with a head bead (cyan) that is linked to four tail beads (crimson) which are again linked to another head bead (Figure 1C right). In this way, both bilayer lipids and bolalipids share the same molecular structure and the same interactions between lipid beads. To decouple the effect of the connected geometry of the bolalipids from that of lipid asymmetry, we assume both head beads of a bolalipid to share the same properties. Bolalipids in archaeal membranes can differ in the number of cyclopentane rings or the branching of the tail and thus in the molecular stiffness (Chong, 2010; Chong et al., 2012). To represent this effect, we added two angular potentials between the second and the fourth and the third and the fifth tail bead with variable strength $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}$. By varying $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}$, we can control the molecular stiffness of the bolalipid molecules and thus model different types of bolalipids. The model is simulated within molecular dynamics implemented in the LAMMPS open-source package (Thompson et al., 2022), and simulations are visualized with OVITO Stukowski, 2010. To include the implicit effect of the surrounding water and to simulate membranes at vanishing tension, we used a Langevin thermostat combined with a barostat that kept the membrane at zero pressure in the $x-y$ plane. Details on the computational model are given in Appendix 1–Section 1.

We first tested the ability of both the bilayer and bolalipid molecules to self-assemble into membranes. To do so, we placed dispersed lipids in a periodic 3D box. For both bilayer and bolalipids, we found that membranes self-assembled over a wide range of lipid interaction parameters (Figure 1D, Figure 1—video 1, Appendix 1–Section 2). We then explored the influence of flexibility in the bolalipid molecule. At small values of the molecular rigidity $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}$, single lipids are flexible and can thus adopt a range of possible conformations. The different conformations can be classified by the angle $\theta$ between the two lipid heads (Figure 1E). Two lipid conformations dominated the conformation distribution in the context of a membrane: the U-shaped conformation with both head beads on the same membrane leaflet ($\theta \approx 0$) and the straight conformation with one head bead in each opposing membrane leaflet ($\theta \approx \pi$) (Figure 1—video 2, Figure 1, Appendix 1–Section 3; for comparison, bilayer membranes in Figure 1—video 3 and rigid bolalipid membranes in Figure 1—video 4). In the self-assembled bolalipid membrane, we marked bolalipids as being in the U-shape conformation if $\theta <\pi /2$ and in the straight conformation otherwise.

To study the phase behaviour of bilayer and bolalipid membranes, we analysed the diffusion of single lipids as a function of the lipid interaction parameters ${ϵ}_{\mathrm{p}}$ and $\omega$ (see Appendix 1–Section 4 and Appendix 1–Section 5). The diffusion constant $D$ exhibited a discontinuity as a function of the temperature, which marks the transition as the membrane moves from the gel phase to the liquid phase. The discontinuity occurred at different values of interaction strength ${ϵ}_{\mathrm{p}}$ and interaction range $\omega$ for bilayer and bolalipid lipids, and it also depended on the values of the molecular stiffness $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}$ (Appendix 1—figure 4). In these simulations, the disintegration of the membrane defined an upper limit to the liquid phase and the transition to the gas phase. Based on this classification, we plotted the phase diagram for bilayer membranes (Figure 1F top left), fully flexible bolalipid membranes ($k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=0$, Figure 1F top right), and rigid bolalipid membranes ($k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=5k_{B}T$, Figure 1F bottom left) as a function of the range of the hydrophobic interaction $\omega$ and the temperature $T_{\mathrm{e}\mathrm{f}\mathrm{f}}$. Membranes made of bilayers (blue) and flexible bolalipid molecules ($k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=0$, orange) behaved similarly under these conditions (Figure 1F bottom right). Strikingly, just as observed in extremophile archaea, as the molecular stiffness of bolalipids increased ($k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=5k_{B}T$, magenta), the liquid region was shifted toward higher temperatures and larger values of the interaction range. This is due to the fact that bolalipid molecules are able to engage in more extensive interactions with partners when in the extended conformation, which helps to stabilize the membrane at higher temperatures.

To explore mechanical properties of bolalipid membranes, we chose $\omega =1.5\sigma$ for the remainder of the work (dashed line in Figure 1F bottom right). Figure 2A shows the phase diagram for bolalipid membranes replotted as a function of temperature and bolalipid rigidity for the chosen interaction range. To be able to compare membranes made of bolalipids of different molecular rigidities, we needed to adjust temperature for each to reach similar fluidities, shown by dashed lines in Figure 2A. We then characterised the conformations of individual bolalipids in flat membranes, measured by the fraction of bolalipids in U-shape conformation $u_{\mathrm{f}}$. We find that for flexible bolalipids ($k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=0$), more than 50% of all lipids are in the U-shape conformation (Figure 2B). This fraction decreases with increasing rigidity of bolalipid molecules and vanishes around $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}\ge 2k_{B}T$, for which almost all bolalipids take up linear conformations.

Figure 2 with 2 supplements see all

Download asset Open asset

This behaviour can be easily captured by considering bolalipids as a two-state system, with straight and U-shaped conformations, as argued before (Appendix 1—figure 1). We assumed that the energy of the straight conformation vanishes $E_{\mathrm{s}}=0$ and the energy of the U-shaped conformation reads $E_{\mathrm{u}}=c_0+c_1{k}_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}$, where, for simplicity, we assumed that the energy linearly depends on $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}$ which is the only relevant energy scale in the system and the two constants $c_0$ and $c_1$, which determine the conformation energy of U-shaped bolalipids. The fraction of bolalipids in U-shape conformation then follows $u_{\mathrm{f}}(k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}})=1/(1+\mathrm{e}\mathrm{x}\mathrm{p}(\beta (c_0+c_1{k}_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}))),$ with $\beta =1/(k_{\mathrm{B}}T)$, as shown by the fit in Figure 2B (dashed grey line, $R^2=0.99$, Appendix 1–Section 6). For the fit, it appears that $c_0<0$, which implies that bolalipids in U-shape conformation are slightly favoured over straight bolalipids at $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=0$ (explored in Appendix 1–Section 6).

It was previously hypothesized that an increasing fraction of bolalipids in straight configuration would increase the membrane rigidity (Vitkova et al., 2020). To determine the membrane rigidity using our model, we assessed the height fluctuation spectrum ($h^2$) of flat membranes in a periodic box (Cooke and Deserno, 2005; May et al., 2007). Interestingly, we found that the original theory of Helfrich, 1973 failed to describe the resulting height fluctuation spectrum. However, the extended theory by Hamm and Kozlov, 2000, which also includes the energetic cost of lipid tilt, successfully captured bolalipid fluctuations (Appendix 1–Section 16). In this case, the resulting height spectrum of the membrane at vanishing membrane tension is given by May et al., 2007

where $q=2\pi n/L$ is the wave number ($n\in \mathbb{Z}$), $L$ is the box size, $\kappa$ is the bending rigidity of the membrane, ${\kappa }_{\theta }$ is the tilt modulus and $l_{\theta }=\sqrt{\kappa /{\kappa }_{\theta }}$ is a characteristic length scale related to tilt, which we call the tilt persistence length. Considering Equation 1, the tilt term is expected to matter if the analysed inverse wave numbers become of the same order as the tilt persistence length $l_{\theta }$. The wave numbers that we analysed correspond to wavelengths that are at least twice the thickness of the membrane ($q<2\pi /(12\sigma )\approx 0.5\sigma$). Thus, the tilt term is expected to noticeably contribute if $l_{\theta }>2\sigma$. For typical bilayer membranes, one finds ${\kappa }_{\theta }=12k_{B}T\mathrm{n}{\mathrm{m}}^{-2}$ (May et al., 2007) and $\kappa =20k_{B}T$ (Deserno, 2015), so $l_{\theta }\sim 1nm\approx 1\sigma$ and, therefore, the tilt term can be neglected as practised before (Cooke and Deserno, 2005). However, when the membrane rigidity increases as we expect for bolalipid membranes, the tilt persistence length $l_{\theta }$ increases and the tilt term in Equation 1 becomes relevant.

By fitting the height spectrum for bolalipid membranes (and bilayer membranes for comparison) (Appendix 1—figure 6), we measured the bending rigidity (Figure 2C) and the tilt persistence length $l_{\theta }$ (Figure 2C inset) as a function of $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}$ (see Appendix 1–Section 7, Figure 2—video 1). With increasing bolalipid molecular rigidity $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}$, the bending rigidity $\kappa$ rose from $8k_{\mathrm{B}}T$ and plateaued at $60k_{\mathrm{B}}T$, showing bolalipid membranes can be very rigid while liquid. Strikingly, the increase in membrane rigidity coincided with U-shaped bolalipids vanishing from the membrane (Figure 2B), which confirmed the hypothesis that straight bolalipids render lipid membranes rigid. At the same time, the tilt persistence length $l_{\theta }$ increases with bolalipid rigidity, starting near zero, crossing the $2\sigma$ threshold at $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=1k_{B}T$ and plateauing at $5\sigma$. Since membrane bending and lipid tilting are two modes of membrane deformations that compete, we conclude that bilayer and flexible bolalipids molecules form flexible membranes that prefer to bend rather than to tilt, while bolalipids in straight configuration form rigid membranes that prefer to tilt rather than to bend. Taken together, bolalipid membranes made of flexible lipid molecules are as flexible as lipid bilayers, adopting U-shaped conformations, where those made of bolalipids in straight configurations are rigid.

Changing membrane curvature alters the area differently in the two membrane leaflets. To adapt to the area difference, we thus expect the fraction of U-shaped bolalipids to change as the membrane curvature changes. Moreover, the results of Figure 2B and C showed that the U-shaped bolalipid fraction and the membrane bending rigidity are correlated. As a result, we predict that the fraction of straight versus U-shaped bolalipids in a membrane will change in response to membrane bending, in a way that makes the bending rigidity of a bolalipid membrane curvature dependent.

To investigate the effect, we measured the bending rigidity of bolalipid membranes in cylindrical shapes as a function of their radii (Harmandaris and Deserno, 2006) (see Figure 2D, Figure 2—video 2, Appendix 1–Section 8). We first noticed that while membrane tubes made of bilayer and flexible bolalipids were stable up to small cylinder radii $R$, almost as small as the membrane thickness itself, we found that membrane made from stiffer bolalipids ($k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=1{\mathit{k}}_{\mathrm{B}}\mathit{T}$) ruptured well before. Strikingly, while for stiffer bolalipids ($k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=1{\mathit{k}}_{\mathrm{B}}\mathit{T}$) the U-shaped bolalipid fraction increased strongly over a short range of the mean membrane curvature ($H=1/(2R)$), we only found a small change in the U-shaped bolalipid fraction of flexible bolalipid membranes ($k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=0$) (Figure 2E). Consequently, since the fraction of U-shaped bolalipid molecules controls the bending rigidity of bolalipid membranes (Figure 2B and C), we found that there is a strong dependency of the bending rigidity $\kappa$ on the membrane mean curvature of stiffer bolalipids ($k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=1{\mathit{k}}_{\mathrm{B}}\mathit{T}$) (Figure 2F).

In an elastic material, the strain modulus holds constant and deformation is reversible. For bolalipid membranes at $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=1{\mathit{k}}_{\mathrm{B}}\mathit{T}$, however, the bending modulus decreases when deformation increases, rendering bolalipid membranes hypoelastic. Fortunately, this dependency on curvature does not invalidate our fluctuation results, where the curvature is small enough that its effect on the bending modulus is negligible (Section 15). In contrast, we did not find that $\kappa$ was curvature-dependent for bilayer or flexible bolalipid membranes ($k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=0$). Taken together, the rigidity of bolalipid membranes is not only controlled by the molecular stiffness of their lipid constituents but also by the emerging geometry of the ensemble of lipids. Since membrane geometry and thus membrane rigidity will change upon membrane deformations, this gives rise to hypoelastic material properties.

Another important material parameter is the Gaussian bending modulus $\overline{\kappa }$, which characterizes the reshaping behaviour of fluid lipid membranes under topological changes (Deserno, 2015). $\overline{\kappa }$ is notoriously difficult to measure since it only becomes detectable when the membrane changes its topological state. Continuum membrane theory, combining shape stability arguments and elasticity, predicts $-\overline{\kappa }/\kappa \in [-0.5,-1]$ (Deserno, 2015), where the former value is expected for incompressible membranes. Indeed, most of the numbers we know for the ratio of the two bending rigidities, many of which were deduced from simulations, lie within this range (Hu et al., 2012). Using the same method as developed by Hu et al., 2012, we determined $\overline{\kappa }$ by measuring the closing efficiency of membrane patches into a sphere (Appendix 1–Section 9). At $T_{\mathrm{e}\mathrm{f}\mathrm{f}}=1.2$, we obtained $\overline{\kappa }=-4.30(22)k_{B}T$ and thus a ratio of $-\overline{\kappa }/\kappa =0.89\pm 0.04$ for bilayer membranes, similar to what has been reported previously (Hu et al., 2012). For flexible bolalipid membranes, we got a slightly smaller value for $\overline{\kappa }=-5.04(37)k_{B}T$. Due to the larger bending modulus, however, flexible bolalipid membranes show a significantly smaller ratio $-\overline{\kappa }/\kappa =0.64\pm 0.04$ ($k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=0$). At larger temperature ($T_{\mathrm{e}\mathrm{f}\mathrm{f}}=1.3$), the ratio can be even smaller $-\overline{\kappa }/\kappa =0.45\pm 0.07$ (Appendix 1−Section 9). The result shows that in addition to the differences in bending rigidities, the ratio of the two bending moduli differs strongly in bilayer and bolalipid membranes.

Archaeal membranes contain varying amounts of bilayer lipids (Oger and Cario, 2013; Tourte et al., 2020). The exact bolalipid/bilayer fraction depends on the growth temperature, with higher levels of bolalipids with increasing temperature (Kim et al., 2019), and higher fraction of cyclopentane rings in the tails (Chong et al., 2012). In order to investigate the effect of different lipid contents on membrane mechanical properties, we wanted to model the archaeal membrane by mixing bilayer lipids into bolalipid membranes. Since in our model, the liquid regions of rigid bolalipid membranes and bilayer membranes do not overlap (Figure 1F bottom right), we picked the temperature $T_{\mathrm{e}\mathrm{f}\mathrm{f}}=1.3$ to minimize fluidity mismatch and we set the molecular rigidity $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=2k_{B}T$ to limit U-shaped bolalipids (Figure 2B). We then measured the diffusion constant $D$ as a function of the fraction of bilayer lipids $f^{\mathrm{b}\mathrm{i}}$ (Appendix 1–Section 5). Interestingly, we found that mixing only 10% bilayer lipids into the bolalipid membrane in gel state is enough to fluidize the membrane (Figure 3A). We then measured the bending rigidity and the tilt persistence length $l_{\theta }$ of flat mixture membranes by analysing the fluctuation spectrum. Both the bending rigidity $\kappa$ (Figure 3B) and the tilt persistence length $l_{\theta }$ decreased non-linearly with the bilayer lipid fraction $f^{\mathrm{b}\mathrm{i}}$ (Figure 3B inset). Taken together, the bolalipid membrane can be substantially softened either through bolalipids acquiring the U-shaped conformation or through addition of bilayer-forming lipids.

Figure 3

Download asset Open asset

To investigate the response of bolalipid membranes to large membrane curvature and topology changes like those induced upon vesicle budding, which regularly occurs in archaea, we simulated membrane wrapping of an adhesive cargo bead (Figure 4A, Figure 4—videos 1–3, Appendix 1–Section 10). Importantly, this provided us with a method to study how lipid organisation is affected by externally imposed membrane curvature and mechanics. We first simulated membrane wrapping at different adhesion energies ${ϵ}_{\mathrm{m}\mathrm{c}}$ between lipid head beads and the cargo until we observed that the membrane wrapped the cargo completely (including membrane fission). Then the minimum adhesion energy, for which a membrane bud completely enveloped the cargo bead, is the onset adhesion energy ${ϵ}_{\mathrm{m}\mathrm{c}}^{\ast }$ (Figure 4A), which we measure as a function of the bolalipid stiffness $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}$ (Figure 4B). For small molecular stiffness $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}$, ${ϵ}_{\mathrm{m}\mathrm{c}}^{\ast }$ first increases linearly with $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}$ before it saturates around $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=3k_{B}T$. We expect that the onset energy is proportional to the membrane bending rigidity ${ϵ}_{\mathrm{m}\mathrm{c}}^{\ast }\propto \kappa$, because the bending energy to wrap a spherical particle is size-invariant (Deserno, 2004; Deserno, 2015). When we increased the bending rigidity, through increasing stiffness of bolalipid molecule $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}$, ${ϵ}_{\mathrm{m}\mathrm{c}}^{\ast }$ increased by a factor of 3 (Figure 4B), suggesting that also $\kappa$ increased by a factor of 3. However, from directly measuring the membrane rigidity from the fluctuation spectrum (Figure 2C), we saw that $\kappa$ increased by a factor of 10. To reconcile these seemingly conflicting observations, we reason that the bending rigidity $\kappa$, similar to Figure 2F, is not constant but softens in the range $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}\in [0,2]k_{\mathrm{B}}T$, upon increasing membrane curvature. This is due to the dynamic change in the ratio between bolalipids in straight and U-shaped conformation. Hence, bolalipid membranes show marked hypoelastic behaviour as they soften during reshaping.

Figure 4 with 4 supplements see all

Download asset Open asset

Through analysing the bolalipid conformations, we found that the membrane was able to curve by increasing the fraction of U-shaped bolalipids in the outer layer of the deformation (Figure 4C, Appendix 1–Section 12). To a lesser degree, we also observed this effect on the inner membrane neck (Appendix 1—figure 14). Remarkably, though, even when lipids are so stiff that there are no more U-shaped bolalipids in the flat mother membrane ($k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=2k_{B}T$), the outer layer of the curved membrane retained a non-negligible fraction of $\approx 10\mathrm{\%}$ U-shaped bolalipids, which in turn decreases $\kappa$ and softens the membrane. However, the softening effect on the membrane, indicated through a constant onset energy for $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}\ge 3k_{B}T$ (Figure 4B), persists even for those very stiff bolalipids. Since for stiff membranes, practically all U-shaped bolalipids are gone (Figure 4C), this suggested that an additional membrane-curving mechanism must be involved.

Looking more closely, at high molecular rigidity ($k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}\ge 2{\mathit{k}}_{\mathrm{B}}\mathit{T}$), we observed the formation of multiple pores on the membrane bud, which we quantified by measuring the time-averaged maximum pore diameter (Figure 4D, Figure 4—video 4, Appendix 1–Section 12). While large pores were not observed in the flat membrane, the diameter of membrane pores around the cargo was found to grow with the increase in bolalipid stiffness. We reasoned that pores form when the energetic cost required to change the bolalipid conformation to release bending stress is larger than the energetic cost of opening a lipid edge surrounding the pore. Hence, for relatively flexible bolalipids, U-shaped bolalipids provide the necessary area difference between the outer and inner layer of the membrane bud and thereby soften the membrane. For stiff bolalipid molecules, however, membrane pores start to form to enable membrane curvature as U-shaped bolalipids become prohibited. Both mechanisms help to explain the discrepancy between ${ϵ}_{\mathrm{m}\mathrm{c}}^{\ast }$ and the bending modulus $\kappa$ obtained by studying membrane fluctuations (Figure 2C).

Having shown that bolalipid membranes can effectively soften also by including some amount of bilayer-forming lipid molecules, we next measured the onset energy ${ϵ}_{\mathrm{m}\mathrm{c}}^{\ast }$ for cargo budding in the membranes formed by mixtures of bolalipids and bilayer-forming lipids, as a function of bilayer lipid head fraction $f_{\mathrm{h}}^{\mathrm{b}\mathrm{i}}$ (Figure 5A, Figure 5—video 1). We found that the onset energy sharply decreases with increasing amount of bilayer forming lipids, and plateaus for 50% bilayer head fraction, where it acquires similar values as in the case of fully-flexible bolalipids (Figure 4B). For small bilayer fractions, U-shaped bolalipids localize almost exclusively on the outer layer of the bud (Figure 5B). As the bilayer fraction increases, there is a steady reduction in the percentage of U-shaped bolalipids in the outer layer in favour of bilayer lipids that take their role in supporting membrane curvature, with U-shaped bolalipids completely vanishing at high fractions of bilayer-forming lipids. The fraction of bilayer lipid head beads initially shows an asymmetry between the preferred outer layer and the penalized inner layer around the cargo (Figure 5C), but eventually approaches $f_{\mathrm{h}}^{\mathrm{b}\mathrm{i}}$ in both layers. Taken together, as the bilayer lipid fraction increases, the role of U-shaped bolalipids in making up the asymmetry between the outer and inner layer is taken over by bilayer lipids. Curiously, the addition of bilayer lipids promotes the formation of U-shaped lipids, both in the flat membrane and even in the inner layer around the bud. This is likely to be explained by the fact that when more bilayer lipids are incorporated, the membrane is less densely packed (Appendix 1—figure 12) and thus U-shaped bolalipids are promoted.

Figure 5 with 1 supplement see all

Download asset Open asset

Importantly, we observed nearly no pores in the membrane bud in mixed membranes, even when we only had a very small fraction of bilayer-forming lipids (Figure 5D). Only as the bilayer fraction increased, did we observe the formation of very small pores in the bud. For the flat mother membrane, however, membrane pores started to form with increasing values of $f^{\mathrm{b}\mathrm{i}}$. They acquired sizes similar to those obtained around the bud in pure bolalipid membranes (Figure 4D). The pore formation in the flat mother membrane is likely promoted because the membrane becomes destabilized by the increasing proportion of bilayer lipids which are close to the gas phase. Taken together, bolalipids can bud porelessly when bilayer-forming lipids, which cause membrane softening, are included.

To study the topological difference between bilayer lipids and bolalipids, we started from a coarse-grained model that was previously developed for lipid bilayer membranes by Cooke and Deserno, 2005. Each lipid is composed of a chain of beads. Bilayer lipids are composed of one head bead connected to a chain of two tail beads. In contrast, bolalipids are composed of two head beads linked to a chain of four tail beads (see Figure 1C). Importantly, bilayer and bolalipids have identical head groups.

Each adjacent pair of beads in a lipid is connected by a finite extensible non-linear elastic bond (FENE). For a given bond length $r$, its potential is the sum of an attractive term and a purely repulsive Lennard-Jones potential that enforces volume exclusion

with $K=30{\mathit{k}}_{\mathrm{B}}\mathit{T}/{\sigma }^2$ and maximum length $R_0=1.5\sigma$ in the first term. We note that in the simulations, our time, distance, and energy units are, respectively, $\tau$, $\sigma$, and ${\mathit{k}}_{\mathrm{B}}\mathit{T}$, and the Boltzmann constant $k_B=1$. Consequently, our unit of mass is given by $1m=1{\mathit{k}}_{\mathrm{B}}\mathit{T}{\tau }^2/{\sigma }^2$. The second term of Equation A1 is given by

where $x=\mathrm{m}\mathrm{i}\mathrm{n}(r,r_{\mathrm{c}})/r_m$, where $r_m$ is where the potential reaches its minimum value $U_{\mathrm{m}}$ and $r_{\mathrm{c}}$ is its cutoff. For bonded beads, we set $U_m=1{\mathit{k}}_{\mathrm{B}}\mathit{T}$, $r_{\mathrm{m}}=r_{\mathrm{c}}=2^{1/6}\sigma$, so that the repulsive part is zero for $r>r_{\mathrm{c}}$.

For non-bonded beads, as a first term, we pick also a purely repulsive form with $r_{\mathrm{m}}=r_{\mathrm{c}}=2^{1/6}\sigma \approx 1.12\sigma$, and parametrise on its strength by scanning the interval $U_{\mathrm{m}}={ϵ}_{\mathrm{p}}\in [0.5,2]{\mathit{k}}_{\mathrm{B}}\mathit{T}$. The inverse of this potential depth is the effective temperature $T_{\mathrm{e}\mathrm{f}\mathrm{f}}=k_{\mathrm{B}}T/{ϵ}_{\mathrm{p}}$ in our system. Following Cooke and Deserno, 2005, we scaled down $r_{\mathrm{m}},r_{\mathrm{c}}$ by $0.95$ for interactions with head beads, to ensure no spontaneous curvature in a membrane leaflet (Cooke and Deserno, 2005; Hu et al., 2012). Accordingly, we note that in our snapshots, each head bead is represented as a sphere of diameter $0.95\sigma$ and each tail bead as a bead of diameter $1\sigma$.

To model lipid rigidity between each consecutive three beads $b_1,b_2,b_3$ in a lipid, we set a harmonic angle potential

with $\theta =\widehat{b_1,b_2,b_3}$. If one of the beads is a head bead (i.e. in a bilayer lipid), we set $K=k_0=5k_{B}T$. In contrast, for the two inner angles in a bolalipid, we set the interval $K=k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}$, where $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}$ is a global constant in the simulation, equal for all bolalipids. In our work, for different simulations, we vary it in the interval $[0,5]{\mathit{k}}_{\mathrm{B}}\mathit{T}$. Therefore, the difference between a bolalipid and two bilayer-forming lipid molecules is the extra bond and two equal angle potentials keeping each three beads connected and aligned, respectively.

To model the (implicit) hydrophobic interaction between lipid tail beads, we add a longer range attractive cosine squared potential

where $\mathrm{c}\mathrm{l}\mathrm{i}\mathrm{p}(x,a,b)=max(min(x,b),a)$ and $r_c$ is set to the cutoff $2^{1/6}\sigma$ of the volume exclusion interaction. The repulsive Lennard-Jones and the attractive potential are combined to make the tail-tail interaction repulsive in the range $[0,r_c]$ and attractive in the range $[r_c,r_c+\omega ]$. We take $\omega$, the attractive range width, as a parameter. As we joined two bilayer lipids to form a bolalipid, in our model bilayer lipids and bolalipids share the same hydrophobic interaction.

We initially verified self-assembly of our lipids into a flat membrane: we placed lipids in a dispersed gas configuration in a periodic 3D cube. We then evolved the system with timestep $\delta t=0.01$ under a Langevin thermostat with relaxation time of $1\tau$ and checked a flat membrane patch eventually formed (Figure 1—video 1).

For the rest of our simulations, we pre-assembled the membrane. This was done by placing lipids far from each other in a hexagonal grid, locking the lipids vertical position and conformation, compressing with volume exclusion until a target area per lipid is reached, followed by energy minimisation and a final check that the lipids formed a horizontal single cluster.

For simulating flat membrane patches of bolalipids, we combined the previously used Langevin thermostat with relaxation time of $1\tau$ with a Nosé–Hoover barostat with relaxation time of $10\tau$. In LAMMPS, this amounts to combining the commands ‘fix langevin’ with ‘fix nph.’ We configured the barostat to set lateral pressure $P_{\mathrm{x}\mathrm{y}}$ to zero by re-scaling the simulation box in the $x$-$y$ plane. We compare this setup to a fixed box length setup, and a NPT ensemble setup, in Appendix 1–Section 17.

For a bolalipid in a flat membrane, we found thresholding on the angle between its first and second half $\theta =\mathrm{\angle }\overrightarrow{b_3,b_1},\overrightarrow{b_4,b_6}$ is sufficient to classify its conformation (Figure 1E, right).

Appendix 1—figure 1

Download asset Open asset

We checked for a few parameters that this angle follows a bimodal distribution (Appendix 1—figure 1) and that indeed bolalipids assume either a straight conformation, one head bead in opposing membrane leaflets, with $\theta \approx \pi$ or a U-shaped conformation, both head beads on the same leaflet, with $\theta \approx 0$. Accordingly, we marked a bolalipid as being in the U-shape if $\theta <\pi /2$ and in the straight conformation otherwise (see Figure 1E for a snapshot).

Appendix 1—figure 2

Download asset Open asset

For pre-assembly of flat membranes containing bolalipids, this implied a choice of which conformation to pick for each lipid, straight or U-shaped. We verified that for the limiting cases of flexible and rigid bolalipids, starting the system with all bolalipids in the straight and U conformation, respectively, resulted in the U-shaped bolalipid fraction $u_{\mathrm{f}}$ equilibrating after simulating at most for ${10}^4\tau$ (Appendix 1—figure 2).

Appendix 1—figure 3

Download asset Open asset

After verifying this, we started our simulations with a fraction of bolalipids derived from assuming a two-state system, with $\mathrm{\Delta }E$ given solely by the central angle potentials (e.g. for U-shaped bolalipids, both angles assume values of $\pi /2$). We checked our system had equilibrated quantitatively by observing that the time series of the fraction of bolalipids in the U-shaped conformation approached a steady state (Appendix 1—figure 3). We still equilibrated for ${10}^4\tau$ before taking measures.

For phase determination (Figure 1F), the membrane was assembled in a horizontal grid of lipids with approximately 25² head beads in each leaflet. For bolalipid membranes, we set the initial fraction of U-shaped bolalipids as we did for simulations depicted in Appendix 1—figure 3. We then ran the simulation, halting it automatically if the box size diverged, at which point we could consider the membrane as being in the gas phase. Otherwise, the simulation ran to completion for ${\mathrm{\Delta }}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}=10\times {10}^3\tau$. This was determined to be sufficient for both the box size $L$ and the relative amounts of lipid species and conformations to equilibrate.

For computing the diffusion constant $D$, we first found the corresponding mean square displacement for each $\mathrm{\Delta }t$ and then fitted $D=⟨{\left|\mathrm{\Delta }x\right|}^2⟩/(4\mathrm{\Delta }t)$. To compute the M.S.D. of a specific lipid in a temporary conformation, we averaged $\left|\mathrm{\Delta }x\right|$ over all possible intervals of duration $\mathrm{\Delta }t$ where the conformation was held. We took care to exclude displacements of lipids floating in the gas phase of the simulation. Appendix 1—figure 4 shows $D$ as a function of temperature $T_{\mathrm{e}\mathrm{f}\mathrm{f}}=k_{\mathrm{B}}T/{ϵ}_{\mathrm{p}}$ for the different membrane types. Appendix 1—figure 4 shows a discontinuity of the diffusion constant $D$ as a function of $T_{\mathrm{e}\mathrm{f}\mathrm{f}}$. We took this discontinuity as marking the transition of the membrane from the gel phase to the liquid phase, and determined it for bilayer and bolalipid membranes for different values of $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}$ (Figure 4). For all membranes and parameters tested, this is equivalent to setting the minimum diffusion constant for a liquid membrane to $D=5\times {10}^{-4}{\sigma }^2/\tau$.

Appendix 1—figure 4

Download asset Open asset

We judged our two-state model fitness in Figure 2B by first making the model linear, expressing it as $\log (1/u_{\mathrm{f}}-1)=x_1{k}_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}+x_0$; then we restricted it to points where $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}\le 2k_{B}T$, or, equivalently, when there were on average $\approx 2$ or more U-shaped bolalipids; with this choices, we obtained $R^2=0.99$.

By splitting the average potential energy between an internal contribution (bonds, angles, and pair interactions between particles in the same molecule) and an external contribution (pair interactions between a molecule and its neighbours), we determined the transition energy from straight to U-shaped bolalipids in detail. We found that this transition lowers the internal potential energy of the bolalipid while increasing its interaction energy. In total, we obtained an energy barrier for the transition of $\mathrm{\Delta }{E}_{\mathrm{s}\to \mathrm{u}}=0.79(1)k_{B}T$. Since the fit indicates, however, that the U-shaped bolalipid conformation is preferred over the straight conformation, we conclude that there must be either an entropic contribution to the free energy or an intermolecular interaction energy favouring U-shaped bolalipids.

In order to measure height fluctuation spectrums (Figures 2C and 3B), we used membranes with 60² head beads per leaflet, with minimum ${\mathrm{\Delta }}_{\mathrm{e}\mathrm{q}}$ set to $20\times {10}^3\tau$; total runtime was ${\mathrm{\Delta }}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}=60\times {10}^3\tau$.

For measuring the bending modulus $\kappa$, we follow the analysis done by Cooke and Deserno, 2005. We simulate a horizontal membrane in a periodic simulation box of dimensions $(l_{\mathrm{x}},l_{\mathrm{y}},l_{\mathrm{z}})$, that is horizontally square with $l_{\mathrm{x}}=l_{\mathrm{y}}=L$. Importantly, we keep membrane tension to a minimum by setting the lateral pressure $P_{\mathrm{x}\mathrm{x}}=P_{\mathrm{y}\mathrm{y}}=0$ via a barostat. As a first equilibration check, we consider the lateral box size $L$ time series. Starting by considering the full series, we measure how much the first and last half differ. For each half, we compute the maximum, minimum and the diameter (max - min). If the relative difference is less than 30%, we consider it equilibrated. Otherwise, we exclude the first frame of the time series and repeat the check.

For each of the $N_{\mathrm{f}}$ remaining frames, we intend to obtain the height field of the membrane $h(x,y)$ within the $x$-$y$ plane. However, our simulations are particle-based, and hence our system is discrete. Therefore, we divide the horizontal plane in a regular $m\times m$ grid, where $m$ is the length of a grid cell. We set $m=40$ by trial and error so that no bins will be empty. We compute $h_b$, the average height of bin $b$. We then apply the 2D discrete Fourier transform in the $x$-$y$ plane, obtaining complex components $h_n$ for $n\in [-m/2,...,0,..,m/2]\times [-m/2,...,0,..,m/2]$, with $h_n=h_{-n}$.

Appendix 1—figure 5

Download asset Open asset

We correct for the binning by multiplying by $\mathrm{sinc}(n_x/m)\mathrm{sinc}(n_y/n)$, which mostly affects the smallest wavelengths (Cooke and Deserno, 2005); we also scale by $⟨L⟩$ obtaining $u=h⟨L⟩$. For each mode number vector $n=(n_x,n_y)$ of $u$, we compute its amplitude squared $|u_n{|}^2$ and phase $\mathrm{\angle }{u}_n$ (examples in Appendix 1—figure 5). For both the amplitude and phase, we then compute the autocorrelation time ${\tau }_c$ in order to get the statistical inefficiency $g=2{\tau }_c+1$. The phase should regularly jump through the endpoints $[0,\pi ]$, so checking it too for stationarity allows to exclude modes whose amplitude is static but which have their phase stuck. By doing this, we find that we can improve the analysis compared to the current state of the art (Ergüder and Deserno, 2021), which only checks the amplitude squared. Between both the amplitude squared and phase components, we take the largest $g$, which then gives us the number of uncorrelated data points as $N_{\mathrm{f}}/g$. We accept a mode as equilibrated if the remaining trajectory contains at least 20 uncorrelated data points. The standard deviation of the mean of $|u_n{|}^2$ must then be scaled by $\sqrt{g}$.

For each spectrum measurement, we performed four simulations with different seeds for the thermal noise. We retained only modes which had equilibrated on all replicas and averaged over modes with the same mode number. We checked that the box size $L$ varied less than 1% between replicas. To not impose an artificial variable window in mode number, we computed the first equilibrated mode for all our simulations and took the maximum equilibrated mode number as a global minimum threshold. We found in our simulations $n\ge 2$. This roughly corresponds to a maximum wavelength cutoff of $32\sigma$ given that our simulation boxes have length $60\sigma$. We used twice the membrane thickness, i.e., $12\sigma$ as minimum wavelength cutoff.

According to continuum membrane theory, at zero membrane tension, the fluctuation spectrum is given by Seifert, 1997

where $\kappa$ is the bending modulus, $L$ is the box size, $q=2\pi n/L$ is the (angular) wave number and $n$ is the mode number (Cooke and Deserno, 2005). While this was adequate for bilayer membranes, we needed to add the tilt term, parametrized by the tilt modulus ${\kappa }_{\theta }$ (see main text Equation 1, May et al., 2007) to get reasonable fits for rigid bolalipid membranes. Then, we obtained

We fitted the data using $N$ measurements of mean and mean standard deviation $(y_i,{\sigma }_i,x_i)$, with $y=⟨{\left|u_n\right|}^2⟩=⟨L{⟩}^2⟨{\left|h_n\right|}^2⟩$ and $x=q=2\pi |n|/⟨L⟩$ to each model $f(x_i)=y_i$. In this case, the different models $f(x)$ are given by Equation A5 and Equation A6, where Equation A5 can be derived from Equation A6 by formally setting ${\kappa }_{\theta }=\mathrm{\infty }$. Typical example fits are shown in Appendix 1—figure 6.

Appendix 1—figure 6

Download asset Open asset

We used the reduced $R^2$ value as an indicator of goodness of fit, defined as $R^2=\sum _i((f(x_i)-y_i)/{\sigma }_i{)}^2/N$. Reasonable values were recognized as $R^2\le 1$. $R^2$ is plotted together with fit results for pure bolalipid membranes and bolalipid/bilayer mixture membranes in Appendix 1—figure 7. In the first row, we show the results of the fit of Equation A5 while in the second row, we used Equation A6. In general, for small values of bolalipid rigidity or large bilayer fraction, the fits without the tilt term were still reasonable. However, as the bolalipid rigidity increased or the bilayer fraction decreased, the fits became worse. Then only fits with the tilt term were reasonable. In addition, we plot the same measures for temperature sweeps for the bilayer at $T_{\mathrm{e}\mathrm{f}\mathrm{f}}=1.2$, the flexible bolalipid and the rigid bolalipid membranes in Appendix 1—figure 8. The results are summarized for all membrane types in Appendix 1—figure 9.

We omitted the error in the main text when less than the unit. Moreover, the low values of tilt modulus will necessarily be accompanied by smaller error bars since the smaller the value, the larger the influence the term will have on the amplitudes.

Appendix 1—figure 7

Download asset Open asset

Appendix 1—figure 8

Download asset Open asset

We note that in our large flat membrane simulations ($L\approx 60\sigma$), flexible bolalipid membranes at $\omega =1.5\sigma$ are stable only at temperatures smaller than $T_{\mathrm{e}\mathrm{f}\mathrm{f}}<1.4$. For $T_{\mathrm{e}\mathrm{f}\mathrm{f}}=1.4$ and presumably above, the membrane folds while shrinking the box until self-contact occurs. This is accompanied by massive oscillations of the box pressure, and can be avoided by halving the timestep. Since these points are near the gas phase and halving the timestep did not significantly change measurements for non-collapsing membranes (Section 17), we simply kept the timestep at 0.01. On the other hand, both the bilayer membranes and rigid bolalipid membranes of same size disassemble by pore formation, followed by simulation box expansion in response to the increased pressure, respectively at $T_{\mathrm{e}\mathrm{f}\mathrm{f}}=1.5$ and 2.3. This explains why in Appendix 1—figure 9 we have more data for bilayer at higher $T_{\mathrm{e}\mathrm{f}\mathrm{f}}$ than for the flexible bolalipids.

We also note that for both the bilayer and the flexible bolalipid membrane, in $T_{\mathrm{e}\mathrm{f}\mathrm{f}}\in [1.2,1.3]$, the tilt modulus varies non-monotonically. It first increases, aligning with the expectation that higher pair potential temperature would reduce order and thus increase the cost of local coordinated tilting, i.e., increasing the tilt modulus ${\kappa }_{\theta }$. However, at the higher $T_{\mathrm{e}\mathrm{f}\mathrm{f}}$ it abnormally decreases.

Comparing the bending moduli of the flexible bolalipid membrane ($\kappa \approx 8k_{B}T$, ${\kappa }_{\theta }\approx 30(10)k_{B}T/{\sigma }^2$) to the bilayer membrane ($\kappa \approx 5k_{B}T$, ${\kappa }_{\theta }\approx 33(20)k_{B}T/{\sigma }^2$, Appendix 1—figure 7B), we find similar results. By systematically investigating the membrane rigidity as a function of temperature, we found that in general flexible bolalipid membranes have a slightly increased rigidity compared to bilayer membranes (Appendix 1—figure 9). We also verified that by increasing $T_{\mathrm{e}\mathrm{f}\mathrm{f}}$, a rigid bolalipid membrane softens in the same manner as bilayer and flexible bolalipid membranes (Appendix 1—figure 9).

Appendix 1—figure 9

Download asset Open asset

Since the method is described elsewhere (Harmandaris and Deserno, 2006), here, we detail specifics. For all measurements, we ran four seeds, each over $20\times {10}^3\tau$, integrating in the NVT ensemble using a Langevin thermostat with damping coefficient $1\tau$ and unit temperature. The membrane was assembled into a cylindrical shape, with the number of heads in each leaflet pre-balanced. The stress tensor was measured at $1\tau$ intervals. The radius was measured from trajectory frames saved every $20\tau$ in the following way. First, we excluded lipids in gas phase by clustering. Then, we computed the centre of mass of the membrane and set it as our origin for the cylinder cross-section $x$-$y$ plane. Then we computed the measured radius as $R=1/⟨1/r⟩$, where the average is over beads and where $r=|x^2+y^2|$ is the distance to the cylinder axial radius. Using $1/r$ instead of $r$ directly compensates to first order for the shell volume $2\pi rdr{l}_z$ dependency on $r$.

After qualitatively checking for quick equilibration of the radius, fraction of U-shaped conformers and the stress tensor components, we dispensed with lengthy equilibration, discarding only the first $1000\tau$ of the trajectory. These observables were then averaged over the rest of the trajectory, with errors determined in a blocking-equivalent manner.

Specifically, in these simulations, we focused on three pairs of parameters, which we name bilayer ($T_{\mathrm{e}\mathrm{f}\mathrm{f}}=1.1$), flexible bolalipids ($k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=0,T_{\mathrm{e}\mathrm{f}\mathrm{f}}=1.2$), and slightly rigid bolalipids at $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=1k_{B}T(T_{eff}=1.4)$. To guarantee correct integration, we lowered the timestep until the ratio $|(P_x+P_y)/(2P_z)|$ was less than 0.1. We also balanced increasing the cylinder length, to improve statistics on the stress tensor measurement, against computational performance and the occurrence of long-term deviations from cylindrical shape at high aspect ratio $R/l_z$. We include the exact parameters and data used in the plots in the main text in Appendix 1—tables 1–3.

Appendix 1—figure 10

Download asset Open asset

Appendix 1—figure 11

Download asset Open asset

We will now, in detail, explain the different behaviour of flexible versus stiffer bolalipid membranes under bending. Membrane bending rigidity in bolalipid membranes decreases dramatically once a small fraction of U-shapes is allowed to form, but then plateaus once this U-shape fraction reaches 20%. In a curved bolalipid membrane, U-shapes must accumulate in the outer leaflet to accommodate for area difference. Together, the bending rigidity non-linear dependence on U-shape fraction and the promotion of U-shapes by curvature explain why, in a membrane made of moderately stiff bolalipids ($k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=1k_{\mathrm{B}}T$), which contain very few U-shapes in the flat state, the bending rigidity of the membrane decreases as curvature increases. While in a membrane made of flexible bolalipid molecules ($k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=0$), where many U-shapes are present in the flat membrane, the bending rigidity does not change with curvature.

Bending rigidity $\kappa$ in flat membranes composed of bolalipids decreases dramatically once a small fraction of U-shapes is allowed to form, but plateaus once more than 20% of U-shaped bolalipids are present. In detail, our data shows that with an increasing bolalipid molecular rigidity $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}$, both the number of U-shaped bolalipids decreases (Figure 2B) and the membrane rigidity $\kappa$ increases (Figure 2C). Thus, the correlation suggests that U-shaped bolalipids soften the membrane, in a non-linear way, where most of the change in membrane bending rigidity happens for U-shaped bolalipid fraction <20% (Appendix 1—figure 11).

Separately, membrane curvature affects the area difference between curved membrane leaflets and thus drives U-shape accumulation. To be specific, a cylindrical membrane with area $A$, mean curvature $H$ and thickness $h$ has the outer leaflet with area $A(1+Hh)$ and the inner leaflet with smaller area $A(1-Hh)$. This can be large, in our simulations up to an area change of $Hh=25\mathrm{\%}$. For pure bolalipid membranes, straight bolalipids occupy the same space in each leaflet. Area difference can then be achieved only by having a different amount of U-shaped bolalipids in each leaflet, which can result in a different U-shape fraction between leaflets and thus asymmetry between leaflets. Appendix 1—figure 10 confirms U-shape head fraction asymmetry that increases with curvature, for both flexible ($k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=0$) and moderately stiff bolalipids ($k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=1k_{B}T$).

Together, these two effects result in membrane softening under curvature for the moderately stiff bolalipids, but constant rigidity for flexible bolalipids (Figure 2F). In detail: for membranes composed of moderately stiff bolalipid molecules ($k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=1{\mathit{k}}_{\mathrm{B}}\mathit{T}$), the U-shape bolalipid head fraction only increases in the outer leaflet, going from 10 to 20% (Appendix 1—figure 10). This is in the high sensitivity region where the bending rigidity is expected to change the most (Appendix 1—figure 11). We hypothesize that the molecular rigidity of a U-shaped bolalipid creates compression on the outer leaflet that stabilizes the membrane curvature and thus causes membrane softening. We suspect that for membranes composed of rigid bolalipids ($k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}>1{\mathit{k}}_{\mathrm{B}}\mathit{T}$), the effect is likely not present due to the absence of U-shape formation even under strong bending.

By contrast, for membranes composed of flexible bolalipids ($k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=0$), the U-shaped bolalipid head fraction changes relatively little from its value for flat membranes (from 50% to respectively 60% and 40% for the outer and inner leaflet, Appendix 1—figure 10). This is in the region where the membrane bending rigidity is expected to respond weakly to U-shape fraction (Appendix 1—figure 11). Additionally, the change is symmetric, so presumably the outer leaflet becomes softer as the inner leaflet becomes stiffer, thus creating opposing effects and only weakly affecting the membrane bending rigidity as a whole. We note that the distinction between the U-shape head fraction that we plot (Appendix 1—figure 10) and U-shape fraction (Appendix 1—figure 11) matters little for this analysis.

In order to determine the Gaussian bending rigidity $\overline{\kappa }$ both of bilayer and flexible bolalipid membranes, we closely followed the procedure as detailed in Hu et al., 2012. As a consequence, we here only describe the general idea and provide the specific parameter values that we used. Details of the simulation protocol can be found in the original publication.

In short, we registered the closing frequency of membrane patches that we initialized as spherical caps. By simulating the closing of membrane patches many times and then repeating the procedure for different degrees of spherical cap opening $x$, we could sample the probability $P(x)$ for precurved caps to close. By fitting $P(x)$ with an analytical expression, the folding probability was related to the dimensionless parameter $\xi$, which itself was related to the cap radius $R$, the membrane line tension $\gamma$, the membrane rigidity $\kappa$, and the Gaussian rigidity $\overline{\kappa }$. Using all (known) parameters, we then determined $\overline{\kappa }$, following Hu et al., 2012

For the folding simulations, we either used 1000 bilayer lipids or 500 flexible bolalipids ($k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=0$) and we simulated at $\omega =1.5\sigma$ and $T_{\mathrm{e}\mathrm{f}\mathrm{f}}=1.2$ or $T_{\mathrm{e}\mathrm{f}\mathrm{f}}=1.3$. Initially, we distributed the lipids along a spherical cap. We allowed the lipids to reach equilibrium while we limited the available space of the lipids between two spherical shells. The full equilibration took $10250\tau$, while the last $100\tau$ were used to introduce different seeds. After the equilibration step, we run the simulation until the precurved membrane patch either fully closed or became flat, determined by the relative shape anisotropy ${\kappa }_{\mathrm{s}}^2$. We used the same cut-offs for ${\kappa }_{\mathrm{s}}^2$ as in Hu et al., 2012. We repeated the simulations for 8 different values of $x$ and 200 seeds each. The radius of the spherical caps $R$ was determined from the cap areas. The cap areas were determined from the averaged number of lipids per cap and the area per lipid. For the line tension, we determined the tension that acts on the open edge of a flat membrane in periodic boundary conditions, following the protocol described in Hu et al., 2012. Using the values for $\xi$, $R$, $\gamma$, and $\kappa$, we then determined $\overline{\kappa }$. All values are summarised in Appendix 1—table 4.

To model a spherical cargo of radius $R_{\mathrm{c}}=8\sigma$ being adsorbed by a membrane, we set up a Lennard-Jones potential between the cargo and lipid beads using Equation A2. We parametrise the strength of the potential $U_m={ϵ}_{\mathrm{m}\mathrm{c}}$, calling ${ϵ}_{\mathrm{m}\mathrm{c}}$ the adsorption energy. With lipid tail beads, the interaction is purely repulsive, with $r_m=r_c=2^{1/6}(8+0.5)\sigma \approx 9.5\sigma$. With lipid head beads, we set up an attractive well by setting $r_c=2^{1/6}(8+0.5)\cdot 1.2\sigma \approx 11.5\sigma$, giving the well a width of $\approx 2\sigma$. This limits this attractive interaction to the head beads of the membrane leaflet closest to the cargo.

For budding simulations (Figures 4 and 5), we first pre-equilibrated for ${10}^4\tau$ a membrane with 60² head beads per leaflet before placing the cargo bead on top of the membrane. To displace any lipids that might be inside the cargo’s volume, we performed a short run in the constant volume ensemble, while scaling from zero to full strength the membrane-cargo interaction. We then ran the simulation at constant pressure until the system was in equilibrium for at least ${\mathrm{\Delta }}_{\mathrm{e}\mathrm{q}}=30\times {10}^3\tau$; the end state should then be either partial or full membrane budding. In practice, total runtime was ${\mathrm{\Delta }}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}\in [60,120]\times {10}^3\tau$. For a significant region of membrane parameters of interest, as we increased ${ϵ}_{\mathrm{m}\mathrm{c}}$ we observed the equilibrium state transition directly from non-budded to a collapsed state with the membrane folded around the cargo in a small simulation box. To impose the existence of a region with budding, we kept all membranes in budding simulations minimally stretched by setting the barostat target lateral pressure to $P_{\mathrm{x}\mathrm{x}}=P_{\mathrm{y}\mathrm{y}}=-0.001{\mathit{k}}_{\mathrm{B}}\mathit{T}/{\sigma }^3$.

Note that it was not possible to measure ${ϵ}_{\mathrm{m}\mathrm{c}}^{\ast }$ for $f^{\mathrm{b}\mathrm{i}}=1$ since budding was always followed by membrane disassembly as the bilayer membrane is close to the gas phase for the used parameters.

Appendix 1—figure 12A shows the area per head bead as a function of $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}$ in the inner and outer layer of the membrane bud and the flat mother membrane.

We can use the area per head as a proxy for tension, by comparing with the values for the flat mother membrane. We note that these measurements excluded membrane within $3\sigma$ of a pore. The plot shows that for bolalipid membranes, while the area per head is non-monotonic in $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}$, it varies by less than $\approx 4\mathrm{\%}$. In contrast, for membrane mixture of bilayer and bolalipids, the area per head as a function of $f^{\mathrm{b}\mathrm{i}}$ monotonically increases (Appendix 1—figure 12), varying by $\approx 10\mathrm{\%}$. By default, both inner and outer leaflets would be relaxed if possible, at the same value of area per head as the flat mother membrane. However, the cargo bead competes with tension and pulls head beads from the outer to the inner layer, thus increasing the area per head for the outer layer and decreasing it for the inner layer. This happens both for relatively flexible bolalipids and mixture membranes with increasing bilayer fraction, which can indeed move head beads between leaflets by, respectively, forming and/or flip-flopping U-shaped bolalipids and flip-flopping bilayer lipids. When this is not possible, the area per head values for inner and outer layer are roughly equal, such as for the mixture membranes with nearly no bilayer lipids at $f^{\mathrm{b}\mathrm{i}}\approx 0.2$. This also happens for bolalipids with $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}>2k_{B}T$, where we note that the area per head for the bud layers are higher than that for the relaxed membrane; this can be understood by considering that at these $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}$ values the bud membrane has pores, whose line tension then competes with tension, stretching the membrane. Lastly, the trend inversion for bolalipids at $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=0.5k_{B}T$ can be understood by considering our choice of $T_{\mathrm{e}\mathrm{f}\mathrm{f}}$ for each $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}$ changes slope at the same point (see Figure 2A).

Appendix 1—figure 12

Download asset Open asset

We developed a pipeline for analysing each frame of membrane simulation trajectories, using the pipeline framework and components from OVITO Stukowski, 2010.

Membrane pore identification

By computing the signed distance to the midplane, we assigned a leaflet to each membrane surface element. We then intersected the membrane surface with the midplane surface, obtaining for each pore a line marking its perimeter. For our purposes, it was sufficient to project each pore perimeter into a least-squares fitted plane and compute the area and perimeter from the projected line. For leaflet area measurements, we considered the two surfaces at equal distance between the membrane surface and the midplane.

For the measurements of pore diameter in Figures 4D and 5D, we first took the ensemble average, i.e., time average with rescaling of std. mean deviation according to autocorrelation, of the total pore area.

We consider a point to be part of a pore if its projection in the midplane membrane is within $3\sigma$ of its surface. We can also assign to a point a leaflet corresponding to its signed distance to the midplane. We apply this location classification to the head bead of each bilayer lipid and each half of a bolalipid. Due to the bolalipid’s symmetry, some combinations of conformation (given by the sign of the angle between their head beads) and location are indistinguishable, while others are transition ephemeral states that for our ensemble measures were not relevant (e.g. U-shapes that have one head bead in each layer). When both head beads are on the same leaflet, the only relevant states are the U-shape, where the heads beads angle is $\le \pi /2$ and the flat state, where the angle is $>\pi /2$. When the head beads are in different leaflets, we get the straight state. Thus, each lipid is assigned a state, some of which have a leaflet; additionally, we consider a lipid to be in a pore if any of it(s) head bead(s) is near a pore.

For simulations where pore formation was significant, namely the simulations with pure bolalipid membranes, we present in the main text conformation measurements that exclude lipids that are in a pore. To clarify that this does not qualitatively change our results, we redid the measurements in Figure 4C, considering exclusively lipids near pores and compared (see Appendix 1—figure 13). Not surprisingly, the only value to change significantly is the fraction of flat bolalipids, which is much larger near pores.

Appendix 1—figure 13

Download asset Open asset

To justify our qualitative observations of bilayer lipids and U-shaped bolalipids, respectively in mixtures and pure bolalipid membranes, enriching at the positive curvature leaflets during budding simulations, we plot the average composition of the axially symmetric profile around the cargo bud. To do this, we chose a trajectory section where the shape of the membrane was roughly static. Then, within a $2000\tau$ time interval, we sampled uniformly 100 frames. For each frame, we took a cylindrical section of the simulation centred on the membrane and axis pointing upwards towards $z$. We assigned to each particle a conformation or species, respectively. Using the cargo bead as the origin for cylindrical coordinates $(r,\theta ,z)$, dropping the angle, and using hexagonal binning on the $r,z$ plane, for each conformation/ species we computed the total number of particles observed. Using these totals, we then thresholded on a minimum number of 40 particles observed per bin and plotted the composition as a fraction of $\left(\mathrm{\#}\text{studied conformation or specie}\right)/\left(\text{total in bin}\right)$ (Appendix 1—figure 14).

Appendix 1—figure 14

Download asset Open asset

In this section, we show that little changes in our work when scaling temperature $T$ instead of our effective temperature $T_{\mathrm{e}\mathrm{f}\mathrm{f}}$. If we compare our phase diagrams when scaling $T_{\mathrm{e}\mathrm{f}\mathrm{f}}$ to scaling simulation temperature $T$, we obtain very similar results (Appendix 1—figure 15).

Naturally, the fraction of U-shapes becomes now a function of temperature, but by keeping to the previously defined $(T_{\mathrm{e}\mathrm{f}\mathrm{f}},k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}})$ line, we maintain the same qualitative disappearance of U-shapes as $k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}$ is increased (Appendix 1—figure 16). The effect on bolalipid fraction and bending rigidity becomes, respectively, less abrupt and strong (Appendix 1—figure 17), but qualitatively we have the same trends as discussed in the main text (Figure 2B and C).

Appendix 1—figure 15

Download asset Open asset

Appendix 1—figure 16

Download asset Open asset

Appendix 1—figure 17

Download asset Open asset

In this section, we show that the maximum curvature imposed on the membrane during fluctuation spectrum measurements is not sufficient to either invalidate the use of the Monge gauge or reach the curvature-dependent regime found for semi-flexible bolalipid membranes (see Figure 2).

We first validate the Monge gauge for describing the membrane surface. We can parameterise a periodic surface without overhangs by a height field $h(x,y)$ relative to a reference plane; we can then perform a Fourier transform in space on this field:

where $n$ are mode number vectors, and to each instantaneous real amplitude $h_n:=2a_n$ corresponds one degree of freedom. $a_0$ is just the average membrane height, kept constant by ensuring the membrane momentum is zero, via constraining the initial thermalisation and the thermostat. The wave vector components are $q_{n,i}=2\pi n_i/L_i$; we used $a_n=a_{-n}$ and for the phase ${\varphi }_n=-{\varphi }_{-n}$.

To apply the Monge gauge means taking several approximations on derivatives of $h$, and thus is only valid if $|\mathrm{\nabla }h|\ll 1$. We expand to obtain:

Since we observed amplitude decreasing with $|n|$, the first two modes, with $|n|=1$, contribute the most to the gradient. We discard the other modes and obtain an upper bound:

How small should be $|\mathrm{\nabla }h|$? From Cooke and Deserno, 2005 Figure 5 (bilayer membrane with $w=1.4$, $ϵ=1$), we find $L=50$, $⟨h^2⟩L^2\approx 500$ for $n=1$ (the point x-axis coordinate matches up with expected $(2\pi /50{)}^2=0.016$). Assuming the average was taken over both modes (with $\overrightarrow{n}=(0,1)$ and $\overrightarrow{n}=(1,0)$) we can write:

While we do not have access to the data to determine the RHS, since the LHS is strictly less or equal to the RHS, we obtain that for Cooke and Deserno the threshold for the norm of the gradient was at 8% or more.

In our sims, we left the modes with $n\le 2$ unequilibrated. Therefore, we cannot use averages directly for these modes, since they might smooth out an otherwise excessive deviation from the Monge Gauge. Instead, we compute $max|\mathrm{\nabla }h|$ for each frame, take the average of each $1000\tau$ block to exclude transient behaviour unlikely to affect membrane properties, and take the maximum of these averages over the full trajectory. We take the maximum over replicas.

In our simulations, we thus obtain, for the bilayer at $T_{\mathrm{e}\mathrm{f}\mathrm{f}}=1.1$, 12%; for flexible bolalipids ($k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=0$), 13%; for stiffer bolalipids ($k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=1k_{B}T$), 6%, and for rigid bolalipids ($k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=5k_{B}T$), 4%.

All of these values are close by 5% at most to the threshold used by Cooke and Deserno, so we are confident to apply the Monge Gauge. Arguably, by forcing the most flexible of these membranes to stay flat (e.g. by fixing the simulation box size to a slightly larger value), we could reduce the deviation at the cost of introducing excessive tension, which would mask the effect of the bending modulus $\kappa$.

Improvements notwithstanding, continuing with the Monge Gauge, we have for the maximum mean curvature the expression:

where in the last equality we restricted ourselves to the two largest wavelength modes.

Since we are interested in the steady state, we repeat the procedure used for obtaining the slope, obtaining, respectively, for each membrane, 0.025 (bilayer), 0.028 (flexible bolalipids), 0.012 (stiffer bolalipids) and 0.01 (rigid bolalipids), in ${\sigma }^{-1}$.

The maximum mean curvature for the stiffer bolalipids ($k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=1$) reaches at most $0.012{\sigma }^{-1}$ in a fluctuation measurement. We can interpolate the data from Figure 2F to obtain that at this mean curvature, the bending modulus is reduced by 7%, compared to its value for a flat membrane. Given this is for the maximum value of membrane mean curvature, we expect a much smaller change on the effective bending rigidity. Therefore, we find our analysis of their spectrums, which assumes a constant bending modulus, valid.

Here, we cover two details regarding our fluctuation simulations and analysis.

Moreover, the relevant spectrum does not depend on how the integration is done

We tested bilayer at $T_{\mathrm{e}\mathrm{f}\mathrm{f}}=1.1$, flexible bolalipids ($T_{\mathrm{e}\mathrm{f}\mathrm{f}}=1.2,k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=0$) and rigid bolalipids ($T_{\mathrm{e}\mathrm{f}\mathrm{f}}=1.8,k_{\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{a}}=5$) using different integrators: zero lateral pressure with overdampened Langevin dynamics, implemented in LAMMPS via fix nph and fix langevin, as done in most of our work; a NPT ensemble, implemented in LAMMPS via fix npt; and Brownian dynamics in a fixed volume, implemented in LAMMPS via fix nve and fix langevin, where the box size is set to the average value in nph + langevin simulations. The resulting measurements of $\kappa$ are at most within 10% of the values reported in the main text (Appendix 1—table 5). The value of $l_{\theta }$ does increase by 30% for this bilayer system, but since it is still below our main text threshold of $2\sigma$ our conclusions do not change.

To start with, we get the same results if we equilibrate first with a barostat, fix the box size, and only then sample the fluctuation spectrum (Appendix 1—figure 19,‘langevin’). It becomes then easy to expect that there is no significant difference for our simulations between using fix nph and fix langevin (Appendix 1—figure 19,‘nph +langevin’) versus fix npt (Appendix 1—figure 19,’npt’). Our barostat simply equilibrates our membrane to near zero tension and then does not significantly contribute to the system. While for bilayer and flexible bolalipid membranes, the first modes (i.e. $\left|\overrightarrow{n}\right|<2$) amplitudes are noticeably above the fitted line for the tilt-including model, these can be reached by including a negative tension term with the negligible value of $-0.05k_{\mathrm{B}}T/{\sigma }^2$ (Appendix 1—figure 19, ‘npt’ and ‘nph +langevin’). We call this value negligible, since it is only roughly 3% of the tension necessary to rupture a membrane in this model (Cooke and Deserno, 2005), and its effect is long vanished at the wavelengths we use for the fits in the main text. We understand this not as an issue with our barostating, but simply as consequence of the $q^{-2}$ dependency on tension. Sufficiently large membranes and small bending modulus $\kappa$ will eventually reveal tension present either due to integration error or the approximation of membrane tension by the projected membrane tension. Nevertheless, to show nothing changes in the spectrum used in the main text, we halved the timestep (setting $\mathrm{\Delta }t=0.005\tau$), which reduces the fitted tension by 3 x and allows a good fit over the entire spectrum with the model that only includes tilt and not tension (Appendix 1—figure 19, ‘timestep halved’). For rigid membranes (e.g. rigid bolalipid membranes), all modes are equilibrated for the shorter duration of $60\times {10}^3\tau$ used in the main text results, and we obtain the same results regardless of the integration scheme used (npt, nph +langevin, or fixed box size with Langevin).

Appendix 1—figure 19

Download asset Open asset

Appendix 1—table 5

	${\chi }^2/N$	$\kappa /{\mathit{k}}_{\mathrm{B}}\mathit{T}$	$\mathrm{\Sigma }/\left(k_{\mathrm{B}}T/{\sigma }^2\right)$	$l_{\theta }/\sigma$
langevin	0.72	8.51±0.12	−0.039±0.0004	0.85±0.05
nph +langevin	0.28	8.71±0.11	−0.050±0.004	0.92±0.04
npt	0.25	8.35±0.09	−0.0492±0.0026	0.86±0.04
nph +langevin (short)	0.17	7.88±0.25	N/A	0.60±0.16
nph +langevin (timestep halved)	0.31	9.13±0.12	N/A	0.97±0.05