Pathways for fast and slow fusion of nanovesicles without membrane rupture

Abstract

Eukaryotic cells continuously remodel their membrane architecture by fusion processes, which are initiated by the adhesion of two membranes and eventually lead to a single membrane with a membrane neck or fusion pore. The fusion of cellular membranes involves membrane proteins but the fusion of biomimetic membranes such as lipid bilayers can be induced by bilayer tension even in the absence of proteins. Tension-induced fusion competes however with membrane rupture, which tends to impair the fusion process. Here, we show by molecular dynamics simulations that nanovesicles enclosed by tensionless and asymmetric bilayers can undergo fusion without rupture and that these fusion processes follow two distinct pathways, a slow and a fast one. Fast fusion starts immediately after an initial point contact between the two vesicles has been established whereas slow fusion occurs only after the vesicles have formed a spatially extended contact area. The two pathways are controlled by the stress asymmetry between the two bilayer leaflets or, equivalently, by the resulting transbilayer torque. Our simulation results have important consequences for the free energy landscapes corresponding to fast and slow fusion of nanovesicles, for experimental studies elucidating these fusion pathways, and for protein-mediated fusion of cellular membranes.

---

1 Introduction

The fusion of bilayer membranes is ubiquitous in all eukaryotic cells, where it is crucial for many cellular processes such as vesicle trafficking between different membrane-bound organelles,¹ exocytosis of neurotransmitters, hormones, and exosomes^2,3 as well as the interactions of enveloped viruses with their host cells.⁴ Membrane fusion starts with the adhesion of two membranes and is completed when the two membranes are connected by a membrane neck or fusion pore. The simplest mechanism for membrane fusion is provided by the mechanical tensions acting within the membranes in their prefusion state.⁵ One major drawback of this fusion mechanism is, however, that tense membranes tend to rupture before they fuse. This competition between membrane fusion and membrane rupture strongly reduces the probability for fusion when induced by bilayer tension.⁶

Here, we show, using coarse-grained molecular dynamics simulations, that even tensionless bilayers of nanovesicles undergo fusion provided their outer bilayer leaflets are stretched whereas their inner leaflets are compressed. We identify two distinct pathways for such fusion processes, a fast and a slow one. Fast fusion starts as soon as the two vesicles come into close contact whereas slow fusion requires the formation of an extended contact area between the two vesicles. Furthermore, both fusion pathways are characterized by vanishing bilayer tension which implies that they avoid membrane rupture.

For simplicity, we focus on the fusion of two identical nanovesicles, which are assembled from the same number of lipid molecules and enclose two spatial regions of equal volume as shown in Fig. 1 and 2 further below. The vesicles then undergo translational and rotational diffusion which leads to their initial contact. Subsequently, the two vesicle bilayers attain a variety of intermediates bilayer states until they fuse and form a single vesicle. Whether they fuse via the fast or slow fusion pathway is primarily controlled by the stress asymmetry between the two bilayer leaflets. This stress asymmetry generates a transbilayer torque⁷ that acts to rotate the bilayers’ cross-sections around their mid-surfaces.

Fig. 1 Time series of cross-sections for fast fusion of two nanovesicles with lipid numbers N_ol = 5600 and N_il = 4500 as well as volume parameter ν = ν₀ = 0.966. The vesicles form an initial contact at t = 0 µs, a hemi-fusion stalk at t = 0.4 µs, and fuse at t = t_f = 0.6 µs. The average fusion time 〈t_f〉 as obtained from nine statistically independent simulations is equal to 〈t_f〉 = 3.1 µs.

Fig. 2 Time series of cross-sections for slow fusion of two vesicles enclosed by tensionless bilayers with N_ol = 5700, N_il = 4400 and ν = 0.966. The vesicles come into contact at t = 0 µs, adhere to each other for 0 < t < 6.2 µs, form a hemifused state for 6.2 µs < t < 20.4 µs, and fuse at t = 23.2 µs, forming a fusion pore or neck close to the contact line. The horizontal arrow at t = 19 µs indicates a hemifusion stalk, the vertical arrows indicate pores across one of the bilayers or across the two adhering bilayers, and the horizontal arrow at t = 23.2 µs points towards the fusion pore. The average fusion time 〈t_f〉 as obtained from nine statistically independent simulations is equal to 〈t_f〉 = 15.8 µs.

Our results are obtained via Dissipative Particle Dynamics (DPD), a coarse-grained molecular dynamics simulation method.^8–10 As explained in the Methods section, DPD enables us to prepare and study nanovesicles with a diameter of about 40 nm and to obtain a reliable statistics for the fusion process. In contrast, the main experimental approach to study the morphology of nanovesicles has been imaging by different variants of electron microscopy.^4,11–13 These experimental methods have one major drawback, however, because they provide only a single snapshot of each vesicle. In contrast, computer simulations with molecular resolution as used here can reveal the molecular dynamics of the lipids and the nanoscopic dynamics of the bilayers, thereby monitoring and elucidating the spatio-temporal remodeling of individual nanovesicles.

Our paper is organized as follows. In the next section, we will describe the assembly of lipid nanovesicles and their characterization in terms of lipid numbers and leaflet tensions. In addition, we will also emphasize that the vesicle volume represents another important control parameter and that the nanovesicles attain a unique reference state with tensionless leaflets. Next, we use time series of simulation snapshots to distinguish fast from slow fusion events. This distinction is further corroborated by a systematic analysis of the different clusters as formed by the headgroups and hydrophobic chains of the two vesicle bilayers. The cluster analysis also identifies specific vesicle states corresponding to hemifusion stalks as well as interbilayer bridges formed by protrusions and splayed lipids. The stochastic character of the fusion process is demonstrated by computing the fusion times of individual simulation trajectories as well as the average fusion time. We also emphasize several important consequences of our simulation results for the free energy landscapes of fast and slow fusion, for possible experimental studies, and for protein-mediated fusion of cellular membranes. At the end, we summarize the results of our paper and give a brief outlook on related topics that can be studied using the same computational approach as described here.

2 Assembly and characterization of nanovesicles

2.1 Assembly of lipids into bilayer membranes

Each nanovesicle is enclosed by a bilayer of lipid molecules, which consists of two leaflets and represents the universal building block of both biological and biomimetic membranes. In the simulations, each nanovesicle is assembled from N_ol lipids in the outer leaflet and N_il lipids in the inner leaflet. To reduce the number of system parameters, all vesicles studied here are assembled from the same total number of lipids as given by N_ol + N_il = 10 100. To study the fusion of two nanovesicles, we focus on two identical nanovesicles, each of which is assembled from the same lipid numbers N_ol and N_il = 10 100 − N_ol. In addition, each vesicle encloses the same number of water molecules and, thus, has the same volume. To start the fusion simulations, we place the two vesicles into close proximity, with their outer leaflets separated by a small distance of the order of the bilayer thickness. The vesicles then undergo translational and rotational diffusion, which typically brings them into close contact after a relatively short simulation time.

2.2 Vesicle volume controlled by osmotic conditions

Each vesicle membrane creates two aqueous compartments, an interior compartment enclosed by the membrane and an exterior compartment outside of the vesicle. Furthermore, lipid bilayers represent semipermeable surfaces, which allow the permeation of water but are impermeable to most ions and even small solute molecules. As a consequence, the volume of the vesicle can be controlled by the concentrations of solutes in the exterior and interior aqueous compartments. In the simulations, the inflation and deflation processes are mimicked by changing the number of water beads, which are initially enclosed by the vesicle. Because the permeation of water across the vesicle membranes is very slow on the time scale of the simulations, the number of enclosed water beads is essentially fixed to its initial value and thus represents another useful control parameter.

2.3 Volume parameter and tensionless bilayers

It will be convenient to measure the vesicle volume in terms of the volume parameter ν, which represents the ratio between the actual number of interior water beads and a fixed reference number of these beads, see Methods section. By changing the vesicle volume and thus the volume parameter ν, we also change the bilayer tension Σ, which vanishes for a certain value ν = ν₀ of the volume parameter. Nonzero values of membrane tension will be given in units of k_BT/d², where k_BT is the thermal free energy and d is the bead diameter, which represents the basic length scale in the simulations. For an inflated vesicle with volume parameter ν > ν₀, the vesicle membrane is subject to a positive bilayer tension, Σ > 0, but experiences a negative bilayer tension, Σ < 0, for a deflated vesicle with ν < ν₀. In the following, we will first focus on vesicles with volume parameter ν = ν₀ and vanishing bilayer tension.

2.4 Leaflet tensions of lipid bilayers

Even though the lipid numbers N_ol and N_il are convenient control parameters for the in silico assembly of nanovesicles, they provide only limited insight into the mechanical and elastic states of the lipid bilayers and their leaflets. Indeed, each leaflet of the bilayer membrane can be stretched or compressed when it is subject to a positive or negative leaflet tension. Furthermore, for a lipid bilayer with a fixed total lipid number, we can identify a unique reference state, in which both leaflets are tensionless.^7,14,15 For total lipid number N_ol + N_il = 10 100 as studied here, the reference state with tensionless leaflets is obtained for and .¹⁴

We denote the outer leaflet tension by Σ_ol and the inner leaflet tension by Σ_il. The sum Σ_il + Σ_ol of the two leaflet tensions is equal to the bilayer tension Σ. The reference state of the bilayer is then characterized by Σ_ol = Σ_il = 0. To avoid membrane rupture, the vesicle membranes are taken to experience only a low bilayer tension. However, even for vanishing bilayer tension, Σ = Σ_il + Σ_ol = 0, 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 and opposite leaflet tension, that is, if the two leaflet tensions satisfy Σ_ol = −Σ_il. To explore the different bilayer states with vanishing bilayer tension, we vary the stress asymmetry ΔΣ = Σ_ol − Σ_il between the two leaflets by changing the lipid numbers N_ol and N_il = 10 100 − N_ol within the outer and inner leaflets for vanishing bilayer tension.

3 Fusion of nanovesicles with tensionless bilayers

Two nanovesicles, which are enclosed by tensionless bilayers, undergo fusion for a sufficiently large and positive stress asymmetry ΔΣ, corresponding to stretched outer and compressed inner leaflets. Intuitively, these fusion processes can be understood as follows. First, within the emerging contact area, the shielding of the hydrophobic bilayer cores from the water is reduced by stretching the two outer leaflets, thereby increasing their areas per lipid. Second, the stretched outer leaflets try to reduce their areas, thereby opening the emerging fusion pore. More precisely, we observe vesicle fusion in all simulations with tensionless bilayers for lipid number N_ol ≲ 5700 and N_il ≳ 4400. Furthermore, each fusion event can be categorized as “fast fusion” or “slow fusion” as we will first demonstrate by visual inspection of the simulation snapshots.

3.1 Fast fusion process

A fast fusion event starts shortly after the vesicles come into local contact for the first time. One example for a fast fusion process is shown in Fig. 1 for lipid number N_ol = 5600 and stress asymmetry ΔΣ = 2.2k_BT/d². In this specific example, the two vesicles fuse at time t = t_f = 0.6 µs. Based on nine replicas, that is, on nine statistically independent simulations, we obtain the average fusion time 〈t_f〉 = 3.1 µs for N_ol = 5600.

3.2 Slow fusion process

Next, we consider two nanovesicles with lipid numbers N_ol = 5700 and N_il = 4400. Compared to the nanovesicles with N_ol = 5600 and N_il = 4500, only 100 lipids corresponding to about 1 percent of the total lipid number have been reshuffled from the inner to the outer leaflet. This small change in the lipid numbers leads to the reduced stress asymmetry ΔΣ = 1.68k_BT/d² and to slow fusion that occurs only after the vesicles have formed a spatially extended contact area. One example for such a slow fusion process is depicted in Fig. 2. In this case, the vesicles fuse at time t = t_f = 23.2 µs. Based on nine replicas, we obtain the average fusion time 〈t_f〉 = 15.8 µs for N_ol = 5700. Thus, by reshuffling about 1 percent of the lipids from the inner to the outer leaflet, the average fusion time exhibits a five-fold increase.

3.3 Fast and slow fusion for other lipid numbers

When we start from the fast fusion process with lipid numbers N_ol = 5600 and N_il = 4500 as displayed in Fig. 1 and further increase the stress asymmetry ΔΣ by reshuffling a certain number of lipids from the outer to the inner leaflet, we further reduce the shielding of the hydrophobic bilayer cores from the water and thus expect to facilitate the fusion process. This expectation is confirmed for lipid numbers N_ol = 5500 and N_il = 4600 as shown in Fig. S1–S3 of the SI. Indeed, by reshuffling a hundred lipids from the outer to the inner leaflet, the average fusion time 〈t_f〉 is reduced from 〈t_f〉 = 3.1 µs for N_ol = 5600 to 〈t_f〉 = 1.4 µs for N_ol = 5500.

On the other hand, when we start from the slow fusion process with lipid numbers N_ol = 5700 and N_il = 4400 as in Fig. 2 and further decrease the stress asymmetry ΔΣ by reshuffling a certain number of lipids from the inner to the outer leaflet, we improve the shielding of the hydrophobic chain region from the water and thus expect to hinder the fusion process. This expectation is confirmed for lipid numbers N_ol = 5900 and N_il = 4200 as shown in Fig. S4 and S5 of the SI. In this case, we reshuffle 200 lipids from the inner to the outer leaflet, which leads to a strong increase of the average fusion time 〈t_f〉 from 〈t_f〉 = 15.8 µs for N_ol = 5700 to 〈t_f〉 > 140 µs for N_ol = 5900. In fact, for N_ol = 5900, we observe only adhesion and no subsequent fusion of the two nanovesicles even for long run times up to t = 140 µs. Therefore, for N_ol = 5900, we obtain only the lower bound 〈t_f〉 > 140 µs for the average fusion time 〈t_f〉.

Our results about the average fusion times for nanovesicles enclosed by tensionless bilayers are summarized in Table S1 of the SI. In addition, we also demonstrate that the distinction between the fast and slow fusion pathways is relatively insensitive to the precise value of the volume parameter ν. Indeed, this distinction also applies to slightly inflated nanovesicles, which exhibit a volume parameter ν = 1 > ν₀ and are subject to a positive bilayer tension, Σ > 0, see Fig. S6–S8 of the SI. Likewise, we also observed slow fusion for a slightly deflated nanovesicle with N_ol = 5700 and volume parameter ν = 0.9 < ν₀, which leads to a negative bilayer tension, Σ < 0, see Fig. S9 of the SI.

Our conclusions about the individual fusion events as obtained by visual inspection of the time series of simulation snapshots can be scrutinized and corroborated in a quantitative manner by counting the total number of clusters, which are formed (i) by the lipid headgroups and (ii) by the lipid chains, as explained in the next subsection.

4 Cluster analysis of simulation trajectories

4.1 Cluster analysis for fast vesicle fusion

Two nanovesicles with outer leaflet numbers N_ol = 5600 undergo fast fusion as displayed by the time series of snapshots in Fig. 1. The corresponding cluster analysis is depicted in Fig. 3, where the number n_H of lipid headgroup clusters as well as the number n_C of lipid chain clusters are shown as a function of time t. In this example, the fusion time t_f = t₂ − t₁ = 0.6 µs, corresponding to the horizontal arrow in Fig. 3.

Fig. 3 Cluster analysis of the fast fusion trajectory in Fig. 1 with N_ol = 5600: Headgroup cluster number n_H (blue) and chain cluster number n_C (red) versus time t. For t < 0 µs, the cluster numbers n_H = 4 and n_C = 2 represent two well-separated vesicles. At time t = t₁ = 0 µs, the two vesicles come into contact. During the time interval 0 µs < t < 0.3 µs, the adhesion of the two vesicles leads to the cluster numbers n_H = 3 and n_C = 2. For 0.4 µs < t < 0.5 µs, we observe hemifusion with n_H = 3 and n_C = 1. Finally, for t ≥ 0.6 µs ≡ t₂, the fused vesicles are characterized by n_H = 2 and n_C = 1. Thus, the fusion time t_f = t₂ − t₁ = 0.6 µs, as shown by the horizontal arrow.

4.2 Hemifusion diaphragms versus hemifusion stalks

Inspection of Fig. 3 reveals that the two vesicles with N_ol = 5600 exhibit an intermediate state with cluster numbers n_H = 3 and n_C = 1 during the time interval 0.4 µs < t < 0.5 µs. These cluster numbers correspond to hemifused states of the two vesicles, in which the outer or proximal leaflets of the two vesicle bilayers have merged whereas their inner or distal leaflets remain separate. In the literature, hemifused states are often described as hemifusion diaphragms.^16–19 Such a diaphragm is provided by a weakly curved bilayer segment which is bounded by a circular line of chain junctions, at which three chain layers meet. In our simulations, we do not observe such hemifusion diaphragms but rather hemifusion stalks as displayed in Fig. 4.

Fig. 4 A hemifusion stalk formed by two adhering vesicle bilayers with N_ol = 5600 lipids at time t = 0.4 µs, see the corresponding snapshot in Fig. 1: Sections across the stalk displaying (a) both the blue headgroups and the red lipid chains of the adjacent bilayers; (b) only the red lipid chains, which form a single chain cluster corresponding to n_C = 1; and (c) only the blue headgroups, which form three headgroup clusters corresponding to n_H = 3. The three headgroup clusters are provided by the merged outer leaflets of the two nanovesicles and by the two inner leaflets of these vesicles.

4.3 Cluster analysis for slow vesicle fusion

Two nanovesicles with outer leaflet numbers N_ol = 5700 undergo slow fusion as illustrated by the snapshots in Fig. 2. The corresponding cluster analysis is depicted in Fig. 5. After the two vesicles come into contact with each other at time t = 0 µs, they adhere and form an extended contact area as in Fig. 2. Such an adhering state is characterized by the cluster numbers n_H = 3 and n_C = 2 as previously shown in Fig. 3 for N_ol = 5600. However, for N_ol = 5700, the adhesion of the two vesicles involves recurrent fluctuations of the chain cluster number n_C between n_C = 2 and n_C = 1.

Fig. 5 Cluster analysis of the slow fusion trajectory in Fig. 2 with N_ol = 5700: Time dependence of headgroup cluster number n_H (blue) and chain cluster number n_C (red). The two vesicles are well separated for t < 0 µs and first come into contact at t = t₁ = 0 µs. For 0 µs < t < 6.2 µs, they adhere to each other with n_H = 3 and n_C = 2. Within this adhering state, the chain cluster number n_C exhibits recurrent fluctuations between n_C = 2 and n_C = 1 until n_C changes permanently from n_C = 2 to n_C = 1 at t = 6.2 µs. At the latter time point, the two vesicles form a hemifused state with n_H = 3 and n_C = 1. This hemifused state persists up to t = 20.4 µs, when the headgroup cluster number n_H changes from n_H = 3 to n_H = 2. A detailed comparison of the cluster analysis with the cross-sections in Fig. 2 reveals that the two vesicles fuse at t = 23.2 µs. The recurrent fluctuations of n_C for t < 6.2 µs are caused by lipid protrusions and lipid flip–flops between the two adhering bilayers, see Fig. 6 and 7, whereas the final changes of n_H for t > 20.4 µs arise from pore formation.

4.4 Interbilayer bridges by protrusions and splayed lipids

The recurrent fluctuations of the chain cluster number n_C between n_C = 2 to n_C = 1 as observed in the adhesion state for the time interval 0 µs < t < 6.2 µs, see Fig. 5, are caused by lipid protrusions and lipid flip–flops between the two adhering bilayers. These lipid protrusions and lipid flip–flops form interbilayer bridges between the adjacent bilayers as shown in Fig. 6 and 7, respectively.

Fig. 6 Protrusions of two highlighted lipids with yellow chains and cyan head groups that form a transient bridge between two adhering bilayers for N_ol = 5700. The protrusion-induced bridge decays quite fast, see the data for the chain cluster number n_C (red) in Fig. 5 with t ≳ 1.1 µs.

Fig. 7 Splayed conformation of a single lipid molecule that forms a transient bridge between the two adhering bilayers for N_ol = 5700. In the splayed conformation, the two chains of the lipid point in opposite directions. The cluster analysis implies that the transient bridge formed by the splayed lipid decays quite fast, see the data for the chain cluster number n_C (red) in Fig. 5 with t ≳ 6.1 µs.

In Fig. 6, we see two highlighted lipids that protrude from the upper and lower bilayer in such a way that the chains of these protruding lipids come into contact. In this way, the chain cluster number n_C is reduced from n_C = 2 to n_C = 1. In Fig. 7, on the other hand, the interbilayer bridge is formed by the splayed conformation of a single lipid molecule, in which the two chains of the lipid point in opposite directions. When the splayed lipid is in contact with the two chain layers of the adhering bilayers, the chain cluster number n_C changes from n_C = 2 to n_C = 1 as in Fig. 5. When the splayed lipid relaxes and returns to the usual lipid conformation with a parallel alignment of its two chains, the cluster number n_C changes back from n_C = 1 to n_C = 2. Splayed lipids, originally denoted as extended lipid conformations,^20,21 have been previously observed in coarse-grained molecular dynamics simulations based on DPD⁶ and on the Martini force field.²²

4.5 Cluster analysis for other lipid numbers

The cluster analysis for fast fusion as observed for nanovesicles with lipid numbers N_ol = 5500 and N_il = 4600 within the outer and inner bilayer leaflets is displayed in Fig. S2 of the SI. Furthermore, the cluster analysis for N_ol = 5900 and N_il = 4200 is shown in Fig. S5 of the SI. As previously mentioned, the two nanovesicles with N_ol = 5900 form an extended contact area for at least 140 µs but are not observed to undergo fusion during this time scale.

5 Statistics of fusion times

In order to fuse, the two vesicles have to adhere and subsequently overcome one or several free energy barriers. As a consequence, the fusion time t_f of an individual fusion event represents a stochastic variable, which can be characterized by a cumulative distribution function P(t_f) and by an average fusion time 〈t_f〉. The cumulative distribution function P(t_f) represents the probability for a fusion event to occur at time t with 0 < t < t_f.

5.1 Fusion time statistics for fast fusion

The fast fusion process of two vesicles with outer lipid number N_ol = 5600 as in Fig. 1 leads to the cumulative distribution function P(t_f) in Fig. 8. This distribution P(t_f) is based on nine simulation trajectories, each of which with a run length of 15 µs. An average over the individual fusion times as measured for the nine trajectories leads to the average fusion time 〈t_f〉 = 3.1 µs. The example in Fig. 1 with t_f = 0.6 µs corresponds to the 1st fusion event of the distribution P(t_f) as displayed in Fig. 8.

Fig. 8 Cumulative probability distribution (black) for the fusion times t_f as obtained from nine fusion simulations of two nanovesicles with N_ol = 5600, N_il = 4500, and ν = 0.966. The blue line represents a fit to the Weibull distribution P (t_f) = 1 − exp (−(t_f/τ)^k) with the fitting parameters τ = 3 µs and k = 1.4. Averaging over the nine discrete fusion times, one obtains the average fusion time 〈t_f〉 = 3.1 µs. The example in Fig. 1 with t_f = 0.6 µs provides the 1st fusion event of the distribution as indicated by the broken vertical line.

5.2 Fusion time statistics for slow fusion

The slow fusion process of two vesicles with outer lipid number N_ol = 5700 as illustrated in Fig. 2 leads to the cumulative distribution function P(t_f) as depicted in Fig. 9. This distribution is obtained from nine simulation trajectories, each of which with a run length of 30 µs. The individual fusion times as measured for the nine trajectories leads to the average fusion time 〈t_f〉 = 15.8 µs. The example in Fig. 2 with t_f = 20.3 µs corresponds to the 7st fusion event of the distribution P(t_f) in Fig. 9, see dashed vertical line in this figure.

Fig. 9 Cumulative probability distribution (black) for the fusion times t_f as obtained from nine fusion simulations of two nanovesicles with N_ol = 5700, N_il = 4400, and ν = 0.966. The blue line represents a fit to the Weibull distribution P (t_f) = 1 − exp (−(t_f/τ)^k) with the fitting parameters τ = 16.2 µs and k = 2.2. By averaging over the nine discrete fusion times, one obtains the average fusion time 〈t_f〉 = 15.8 µs. The example in Fig. 2 with t_f = 20.3 µs provides the 7th fusion event of the distribution as indicated by the broken vertical line.

5.3 Fusion time statistics for other lipid numbers

For lipid numbers N_ol = 5500 and N_il = 4600 as well as volume parameter ν = ν₀ = 0.966, the two nanovesicles undergo a fast fusion process, see Fig. S1 and S2 in the SI. The corresponding fusion time statistics is displayed in Fig. S3 of the SI.

6 Important consequences of simulation results

6.1 Free energy landscapes

The fusion process observed in our simulations represents a topological transformation that is not coupled to any chemical reaction such as nucleotide hydrolysis. As a consequence, this process is “downhill” or exergonic in the corresponding free energy landscape, which has two local free energy minima corresponding to (i) the initial two-vesicle state before fusion (BF) and (ii) the final one-vesicle state after fusion (AF). In order to proceed without any chemomechanical coupling, the free energy of the two-vesicle BF state must lie above the free energy of the one-vesicle AF state. Furthermore, the rate, with which the fusion process proceeds, is governed by the free energy barrier that separates the BF state from the AF state. Indeed, the fusion pathways observed here imply that fast fusion has to overcome only a relatively low barrier whereas slow fusion encounters a relatively high barrier. Our simulations also demonstrate that the barrier height can be controlled by the stress asymmetry between the two bilayer leaflets, with a larger stress asymmetry generating a lower free energy barrier.

It is also instructive to compare the curvature energies of the initial BF states and the final AF states. These curvature energies consist of two terms, the bending energies and the Gaussian curvature energies.^23,24 Inspection of Fig. 1 for fast fusion implies that the bending energies of these two states are very similar because the membrane neck in the third snapshot of Fig. 1 does not make any significant contribution to the bending energy. Therefore, the curvature energies of these two states are determined by the Gaussian curvature energies which are proportional to the Gaussian curvature (or saddle-splay) modulus κ_G. More precisely, the two-vesicle BF state has the Gaussian curvature energy 8πκ_G whereas the one-vesicle AF state has the Gaussian curvature energy 4πκ_G. The “downhill” fusion process from BF to AF then implies that the free energy change is given by −4πκ_G, which is negative for positive κ_G. For slow fusion as in Fig. 2, the free energy of the vesicle system involves adhesion free energies in addition to the curvature energies. As far as the curvature energies are concerned, the transformation of the BF state at time t = 0 μs to the AF state at t = 23.2 μs is again dominated by the change −4πκ_G of the Gaussian curvature energies as follows from a comparison of the corresponding snapshots in Fig. 2. Furthermore, the adhesion free energy of the AF state acts to further decrease the free energy of the one-vesicle AF state and thus to stabilize the latter state compared to the two-vesicle BF state.

For symmetric planar bilayers, the Gaussian curvature modulus κ_G can be obtained from the second moment of the stress profile.²⁵ Using previous results for the stress profile of such a bilayer,⁷ we have confirmed that κ_G is indeed positive and has the numerical value κ_G = (7 ± 1) k_BT for the molecular lipid model used here. A positive value κ_G ≃ 20k_BT has also been reported for the Gaussian curvature modulus of egg lecithin bilayers at room temperature, based on the experimental observationave that these bilayers form a lattice of passages.^25,26 The same positive value κ_G ≃ 20 k_BT was also obtained from the theoretical analysis²⁷ of electron microscopy images for smectic lyotropic phases as formed by egg lecithin bilayers.²⁸ On the other hand, giant vesicles enclosed by lipid bilayers that contain two phospholipids and cholesterol have been observed to divide into two daughter vesicles without any chemomechanical coupling, which implies that these lipid bilayers have a negative Gaussian curvature modulus.²⁹

6.2 Experimental studies of nanovesicle fusion

Nanovesicles as investigated here by molecular dynamics simulations are also accessible to experimental studies. Indeed, both synthetic and cellular nanovesicles have been experimentally studied for a long time. Synthetic nanovesicles are assembled from lipid molecules, using a variety of preparation methods, which produce a wide range of vesicle sizes.³⁰ These methods include extrusion of lipid dispersions through filters with a certain pore size^31,32 and microfluidic mixing.³³ In vivo, even smaller nanovesicles are frequently observed such as synaptic vesicles with a diameter between 20 and 60 nm^34,35 as well as exosomes with a diameter between 40 and 100 nm.^36,37

Furthermore, several experimental methods have been developed to change the molecular densities within the two leaflets of the nanovesicles. These methods include: asymmetric insertion of amphiphilic molecules with a slow transbilayer flip–flop rate;³⁸ fusing the vesicles with micelles that carry ligand molecules;³⁹ cyclodextrin-mediated exchange of lipids between two vesicle populations;⁴⁰ flip–flops of an anionic phospholipid from the outer to the inner leaflet driven by Ca²⁺ ions in the exterior solution;^41,42 and enzymatic conversion of phospholipid headgroups in the outer leaflet.^43,44 All of these protocols change the molecular densities in the bilayer leaflets and, thus, can be used to vary the leaflet tensions.

6.3 Local stress asymmetry generated by membrane proteins

Fusion of cellular membranes also proceeds without membrane rupture but typically involves the formation of protein complexes anchored to the membranes in their pre-fusion state. Relatively simple examples are provided by viral fusion proteins such as hemagglutinin molecules⁴⁵ that are anchored to the membrane envelope of the virus. The viral fusion protein inserts into the acceptormembrane of the host cell, a process that does not seem to require additional proteins within the acceptor membrane. Another protein-mediated fusion mechanism is used in homotypic fusion of intracellular membranes that enclose the endoplasmic reticulum or belong to the outer membranes of mitochondria. In this case, the apposing membranes are pulled together by the dimerization of identical membrane proteins, which are anchored in both membranes to be fused. These membrane proteins have a GTPase domain that hydrolyzes GTP. Examples are proteins from the atlastin family, which drive GTP-dependent fusion of the endoplasmic reticulum membranes in multicellular animals.^46–49 A third type of multiprotein complex is used along the outward secretory or exocytic pathway in eukaryotic cells. This pathway, which involves the heterotypic fusion of different organelles and vesicles, is dominated by proteins from the SNARE family.^50–52 For SNARE-mediated fusion, each of the fusing membranes contributes a different set of membrane proteins to the multiprotein complex that mediates their adhesion. Based on our simulation results for nanovesicles, we propose that these membrane proteins generate a local stress asymmetry together with an associated torque that leads to a stretching of the proximal leaflets and to a compression of the distal leaflets of the two membranes in their prefusion state.

7 Summary and outlook

In this paper, we use dissipative particle dynamics (DPD), a coarse-grained molecular dynamics method, to show that two lipid vesicles can undergo tension-induced fusion without membrane rupture. All fusion events as observed in our simulations can be classified as “fast fusion” or “slow fusion”. We first analyze these fusion events via time series of simulation snapshots, see Fig. 1 and 2, using the lipid numbers N_ol and N_il assembled in the outer and inner leaflets of the two vesicle membranes as convenient control parameters. More precisely, the vesicles are assembled from N_ol lipids in the outer leaflets and N_il = 10 100 − N_ol lipids in the inner leaflets of the vesicle bilayers.

In the main text, we study tensionless bilayers with volume parameter ν = ν₀. Fast fusion is observed for N_ol = 5600, as illustrated in Fig. 1 whereas slow fusion dominates for N_ol = 5700 as depicted in Fig. 2. For N_ol = 5600, eight out of nine simulation trajectories lead to fast fusion with an average fusion time 〈t_f〉 = 3.1 µs (Fig. 8). In contrast, for N_ol = 5700, that is, after reshuflling only 100 lipids, which represent less than one percent of the total lipid number, from the inner to the outer leaflet, eight out of nine simulation trajectories lead to slow fusion with an average fusion time 〈t_f〉 = 15.8 µs (Fig. 9). Thus, reshuffling less than one percent of the lipids from the inner to the outer leaflet leads to a five-fold increase in the average fusion time, corresponding to an increase in the free energy barrier of about 1.6 k_BT between the BF and AF states, assuming an Arrhenius-type kinetics. Likewise, fast fusion is also observed for N_ol = 5500 (Fig. S1–S3 in the SI) whereas for N_ol = 5900, even rather long simulation trajectories with run times above 140 µs do not lead to vesicle fusion (Fig. S4 and S5 in the SI), which implies that the individual fusion times exceed 140 µs for N_ol = 5900.

To corroborate the results as obtained from the time series of simulation snapshots, we study the individual simulation trajectories by cluster analysis, that is, by counting the numbers n_H and n_C of hydrophilic headgroup and hydrophobic chain clusters as a function of time. This cluster analysis allows us to identify transient states of the interacting vesicles as provided by hemifusion stalks (Fig. 4), lipid protrusions (Fig. 6) and splayed lipids (Fig. 7). Both lipid protrusions and splayed lipids reduce the number n_C of chain clusters, see the cluster analysis in Fig. 5. Finally, we looked at the statistics of membrane fusion and computed the cumulative probability distribution P (t_f) for the fusion time t_f as displayed in Fig. 8 and 9 for different bilayer asymmetries.

We also emphasized some important consequences of our simulation results as discussed in the previous subsection. First, the free energy landscape for this fusion process is characterized by a positive value of the Gaussian curvature modulus. Second, we describe how the stress asymmetry between the two leaflets of the nanovesicles could be changed experimentally, using a variety of preparation methods that have been discussed in the literature. Third, we propose that protein-mediated fusion of cellular membranes without membrane rupture is accomplished by local stress asymmetries within the bilayer segments close to the proteins, thereby stretching the proximal leaflets and compressing the distal leaflets of the bilayer segments to be fused.

For the slow fusion pathway, we observed membrane pores forming and closing spontaneously, as indicated by the vertical arrows in Fig. 2 and Fig. S8 and S9. We emphasize that these transient pores are distinct from irreversible membrane rupture as observed in previous simulation studies.^5,6 In fact, both the fast and slow fusion processes as described here occur without rupture of the vesicle membrane.

It is worth noting that during vesicle fusion, lipids can exchange between the inner and outer leaflets of the two bilayers, causing the lipid number asymmetry to be altered and relaxed to some extent. This lipid exchange is facilitated by the disordered bilayer structure within the contact area of the two vesicles. As a consequence, the initial asymmetry of the outer and inner leaflets consisting of N_ol and N_il lipids is not preserved in the bilayer of the final vesicle formed by the fusion processes. In particular, the bilayer of the final vesicle is not necessarily tensionless. Moreover, the lipid exchange during vesicle fusion is a stochastic process, and the tensions of the leaflets in the final vesicle depend on the particular realization of this process.

Our computational approach can be extended to a variety of membrane systems. One may study, for instance, how fast and slow fusion are affected by the geometry of the vesicles, that is, by vesicle size and membrane curvature. Using coarse-grained molecular dynamics simulations, such a dependence on the vesicle geometry has been previously observed for lipid flip–flops between the two bilayer leaflets as well as for structural instabilities of the lipid bilayers, revealing a strong increase of the flip–flop rates with increasing stress asymmetries.^15,53 Very recently, the activation energies for lipid flip–flops have also been studied using a model for glassy dynamics.⁵⁴

Likewise, it will be interesting to consider the limit, in which one of the vesicles becomes a planar bilayer with vanishing curvature. Irrespective of its size and curvature, each vesicle will attain a unique reference state with tensionless leaflets as well as tensionless bilayer states with one stretched and one compressed leaflet. On the other hand, when the vesicle bilayers are subject to different bilayer tensions, they will exhibit shape fluctuations of different amplitudes, which are expected to affect their fusion statistics as well.

For simplicity, our study of vesicle fusion focuses on bilayers with a single lipid component. For bilayer membranes with two or more lipid components, one can study the fusion of nanovesicles with different lipid compositions. Furthermore, multicomponent vesicle membranes can undergo lipid phase separation and form intramembrane domains of different compositions. The vesicles will then adhere and fuse via specific domains with a distinct lipid composition. Finally, our study is based on the DPD simulation technique, which represents a coarse-grained molecular dynamics method with particularly simple force fields. However, the distinction between fast and slow fusion should also be accessible to other molecular dynamics simulations with more complex force fields that are derived from all-atom force fields in a more systematic manner. Such more elaborate force fields can be used to study how fast and slow fusion processes depend on the chemical structure of the lipids.

8 Methods

8.1 Coarse-grained model for lipid molecules

Lipid molecules are modeled as in our previous studies^7,14,53: Each lipid has a head group, consisting of three hydrophilic beads, as well as two chains, each of which consists of six hydrophobic beads. All beads have the same diameter d, which is about 0.8 nm and provides the basic length scale for the simulations.

8.2 Assembly of nanovesicles

The software package Packmol⁵⁵ is used to assemble two spherical vesicles, each with N_il lipids in the inner leaflet and N_ol lipids in the outer leaflet. In all simulations described here, the total number of lipids is always taken to be N_il + N_ol = 10 100 for each vesicle. In this way, we reduce the lipid number space from two to one dimension because the outer leaflet number N_ol determines the inner leaflet number. Each vesicle has the outer radius 24 d, corresponding to a diameter of about 40 nm. The number of water beads inside each vesicle is N_W ≡ νN^isp_W with N^isp_W ≡ 90 400 and the volume parameter ν < 1. The simulation box is a cuboid of side lengths L_x = L_y = 86 d and L_z = 130 d. Outside the vesicles, the box contains about 2.5 million water beads. Initially, the Cartesian coordinates of the centers of the vesicles are (0,0,25d) and (0,0,−25d). Thus the initial distance between the centers of the vesicles is 50 d. Moltemplate⁵⁶ is used to convert Packmol output files into the LAMMPS format.

8.3 Dissipative particle dynamics

The fusion of two nanovesicles is studied by dissipative particle dynamics (DPD),^8–10 a coarse-grained simulation method of molecular dynamics. Compared to all-atom molecular dynamics simulations, the DPD method has several advantages. First, one can study the behavior of relatively large bilayer segments, which are not accessible to all-atom simulations. Second, the molecular models used in DPD are built up from a small number of different molecular groups or beads and, thus, involve only a small number of parameters for the interaction forces between the beads. As a consequence, one is able to explore large regions of the parameter space and to study the relative importance of these parameters in a systematic manner. Third, DPD simulations explicitly include water, the universal solvent for all biomolecules, and reproduce the correct hydrodynamics, because all forces used in DPD conserve momentum. The hydrodynamic interactions between different membrane segments of the nanovesicles are important because their remodeling involves transient hydrodynamic flows, which affect the time evolution of the systems. Fourth, membrane fusion is slowed down by free energy barriers. The time to cross such a barrier is stochastic in nature. In such a situation, one should sample a sufficiently large number of fusion events in order to obtain a reliable statistics, which would be prohibitively expensive for all-atom molecular dynamics simulations.

8.4 Vesicle simulations

The simulations are performed using LAMMPS version 2 August 2023⁵⁷ with the DPD interaction parameters introduced in our earlier studies on lipid bilayers and vesicles.^7,14 For each vesicle, equilibration simulations of one microsecond are performed using the Berendsen barostat with the damping parameter of 1 ns to reach the standard DPD pressure of 23.7 k_BT/d. Next, production simulations are performed with fixed side lengths of the simulation box. The time step of 10 ps is used in both the equilibration and production runs. Simulations in which no contact between the vesicles was observed in 15 µs were restarted with a different seed for the random number generator.

8.5 Cluster analysis of headgroup and chain clusters

To follow the topological changes of the two vesicles, we count the number of different clusters, which are formed by the lipid headgroups and by the lipid chains of the two vesicles. For a given configuration of the two vesicles, say vesicle I and vesicle II, two headgroup beads with indices i and j are assigned to the same cluster if their distance d_ij < 1.5d where d is the bead diameter of the coarse-grained molecular model as previously mentioned. Likewise, two lipid chains belong to the same chain cluster if the two chains involve two chain beads k and l with distance d_kl < d. A similar cluster analysis was previously used for polymers interacting with membranes.⁵⁸

As long as the two vesicles are well separated from each other, the cluster analysis leads to n_H = 4 headgroup clusters and n_C = 2 chain clusters: two headgroup clusters formed by the inner and outer leaflets of vesicle I, two headgroup clusters provided by the inner and outer leaflets of vesicle II, as well as two chain clusters corresponding to the hydrophobic chains the two bilayer membranes. When the vesicles adhere to each other, the cluster number n_H of lipid headgroups is reduced from n_H = 4 to n_H = 3, corresponding to the inner leaflet of vesicle I, the inner leaflet of vesicle II, and the combined outer leaflets of vesicles I and II, whereas the cluster number n_C of lipid chains remains at n_C = 2. After the fusion of the two vesicles, we obtain a single fused vesicle with n_H = 2 clusters of headgroups, corresponding to the inner and outer leaflet of the fused vesicle, and a single chain cluster with n_C = 1.

8.6 Computing fusion times via cluster analysis

Keeping track of how the cluster numbers n_H and n_C change with time t helps to identify characteristic events for a given trajectory, as illustrated by Fig. 3 and 5 for N_ol = 5600 and 5700, respectively. If the vesicles get into contact at time t = t₁, the headgroup cluster number n_H is reduced from n_H = 4 to n_H = 3 whereas the lipid chain cluster number n_C remains at n_C = 2. When the vesicles fuse at time t = t₂, the cluster numbers change to n_H = 2 and n_C = 1. Therefore, we define the fusion time t_f of an individual simulation trajectory by t_f ≡ t₂ − t₁. In all cases, the results of the cluster analysis have been confirmed by visual inspection of the simulation snapshots.

Conflicts of interest

There are no conflicts to declare.

Data availability

The compuer simulations were performed using LAMMPS software. The code for LAMMPS can be found at https://www.lammps.org/download.html with a detailed user guide at https://doi.org/10.1016/j.cpc.2021.108171. The version of the code employed for this study was version 2 August 2023.

Supplementary information (SI) contains Table S1 and Fig. S1–S9. See DOI: https://doi.org/10.1039/d6sm00288a.

Acknowledgements

R. L. thanks Vahid Satarifard for useful discussions. This research was supported by the National Science Centre of Poland via grant number 2021/40/Q/NZ1/00017 as well as by the Max Planck Society and the Dieter Schwarz foundation via the Max Planck School “Matter to Life”. The numerical simulations were carried out using the supercomputer resources at the Centre of Informatics - Tricity Academic Supercomputer and networK (CI TASK) in Gdansk, Poland. Open Access funding is provided by the Max Planck Society.

References

1. J. S. Bonifacino and B. S. Glick, Cell, 2004, 116, 153–166.
2. B. S. Joshi, M. A. de Beer, B. N. G. Giepmans and I. S. Zuhorn, ACS Nano, 2020, 14, 4444–4455.
3. L. Cui, H. Li, Y. Xi, Q. Hu, H. Liu, J. Fan, Y. Xiang, X. Zhang, W. Shui and Y. Lai, Mol. Biomed., 2022, 3, s43556.
4. S. M. Kephart, N. Hom and K. K. Lee, Trends Biochem. Sci., 2024, 49, 916–931.
5. J. Shillcock and R. Lipowsky, Nat. Mater., 2005, 4, 225–228.
6. A. Grafmüller, J. C. Shillcock and R. Lipowsky, Biophys. J., 2009, 96, 2658–2675.
7. B. Różycki and R. Lipowsky, J. Chem. Phys., 2015, 142, 054101.
8. P. Hoogebrugge and J. Koelman, Europhys. Lett., 1992, 19, 155–160.
9. R. Groot and P. Warren, J. Chem. Phys., 1997, 107, 4423–4435.
10. J. C. Shillcock and R. Lipowsky, J. Chem. Phys., 2002, 117, 5048–5061.
11. N. G. Kotla, B. Chandrasekar, P. Rooney, G. Sivaraman, A. Larrañaga, K. V. Krishna, A. Pandit and Y. Rochev, ACS Biomater. Sci. Eng., 2017, 3, 1262–1272.
12. A. L. Robson, P. C. Dastoor, J. Flynn, W. Palmer, A. Martin, D. W. Smith, A. Woldu and S. Hua, Front. Pharmacol., 2018, 9, 10.3389–10.3397.
13. C. Chen, C. Chen, Y. Li, R. Gu and X. Yan, Fundam. Res., 2023, 3, 488–504.
14. R. Ghosh, V. Satarifard, A. Grafmüller and R. Lipowsky, Nano Lett., 2019, 19, 7703–7711.
15. A. Sreekumari and R. Lipowsky, Soft Matter, 2022, 18, 6066–6078.
16. S. Hui, T. Stewart, L. Boni and P. Yeagle, Science, 1981, 212, 921–923.
17. S. Kozlovsky, L. V. Chernomordik and M. M. Kozlov, Biophys. J., 2002, 83, 2634–2651.
18. L. K. Tamm, J. Crane and V. Kiessling, Curr. Opin. Struct. Biol., 2003, 13, 453–466.
19. J. M. Gardner and C. F. Abrams, Biochim. Biophys. Acta, Biomembr., 2018, 1860, 1452–1459.
20. P. K. J. Kinnunen and J. M. Holopainen, Biosci. Rep., 2000, 20, 465–482.
21. H. A. Scheidt, K. Kolocaj, D. B. Konrad, J. A. Frank, D. Trauner, D. Langosch and D. Huster, Biochim. Biophys. Acta, Biomembr., 2020, 1862, 183438.
22. Y. Smirnova, S. Marrink, R. Lipowsky and V. Knecht, JACS, 2010, 132, 6710–6718.
23. W. Helfrich, Z. Naturforsch., 1973, 28c, 693–703.
24. R. Lipowsky, The Giant Vesicle Book, Taylor & Francis, 2019, ch. 5, pp. 73–168.
25. W. Helfrich, Physics of Defects, North-Holland Publishing Company, Amsterdam, 1981, pp. 715–755.
26. W. Harbich, R. Servuss and W. Helfrich, Z. Naturforsch., 1978, 33a, 1013–1017.
27. A. G. Petrov, M. D. Mitov and A. Derzhanski, Phys. Lett., 1978, 65A, 374–376.
28. M. Kleman, C. Williams, M. Costello and T. Gulik-Krzywicki, Philos. Mag., 1977, 35, 33–56.
29. J. Steinkühler, R. L. Knorr, Z. Zhao, T. Bhatia, S. Bartelt, S. Wegner, R. Dimova and R. Lipowsky, Nat. Commun., 2020, 11, 905.
30. B. S. Pattni, V. V. Chupin and V. P. Torchilin, Chem. Rev., 2015, 115, 10938–10966.
31. M. J. Hope, M. B. Bally, G. Webb and P. R. Cullis, Biochim. Biophys. Acta, 1985, 812, 55–65.
32. R. C. MacDonald, R. I. MacDonald, B. P. M. Menco, K. Takeshita, N. K. Subbarao and L.-R. Hu, Biochim. Biophys. Acta, 1991, 1061, 297–303.
33. A. Jahn, S. M. Stavis, J. S. Hong, W. N. Vreeland, D. L. DeVoe and M. Gaitan, ACS Nano, 2010, 4, 2077–2087.
34. L. Qu, Y. Akbergenova, Y. Hu and T. Schikorski, J. Comp. Neurol., 2009, 514, 343–352.
35. B. R. Hopiavuori, L. D. Bennett, R. S. Brush, M. J. V. Hook, W. B. Thoreson and R. E. Anderson, OCL: Oilseeds Fats, Crops Lipids, 2016, 23, D113.
36. R. M. Johnstone, M. Adam, J. R. Hammond, L. Orr and C. Turbide, J. Biol. Chem., 1987, 262, 9412–9420.
37. M. Simons and G. Raposo, Curr. Opin. Cell Biol., 2009, 21, 575–581.
38. L. Hua, S. Akcesme, K. Müller and H. Heerklotz, J. Phys. Chem. B, 2025, 129, 1260–1273.
39. T. Allen, P. Sapra and E. Moase, Cell. Mol. Biol. Lett., 2002, 7, 217–219.
40. H.-T. Cheng, Megha and E. London, J. Biol. Chem., 2009, 284, 6079–6092.
41. 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.
42. H.-Y. Guo, H.-Y. Sun, G. Deng, J. Xu, F.-G. Wu and Z.-W. Yu, Langmuir, 2020, 36, 12684–12691.
43. 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.
44. R. Takaoka, H. Kurosaki, H. Nakao, K. Ikeda and M. Nakano, Biochim. Biophys. Acta, Biomembr., 2018, 1860, 245–249.
45. S. Boonstra, J. S. Blijleven, W. H. Roos, P. R. Onck, E. van der Giessen and A. M. van Oijen, Annu. Rev. Biophys., 2018, 47, 153–173.
46. G. Orso, D. Pendin, S. Liu, J. Tosetto, T. J. Moss, J. E. Faust, M. Micaroni, A. Egorova, A. Martinuzzi, J. A. McNew and A. Daga, Nature, 2009, 460, 978–983.
47. J. Hu, Y. Shibata, P.-P. Zhu, C. Voss, N. Rismanchi, W. A. Prinz, T. A. Rapoport and C. Blackstone, Cell, 2009, 138, 549–561.
48. N. Wang and T. A. Rapoport, J. Cell Sci., 2019, 132, jcs227611.
49. R. Lipowsky, S. Pramanik, A. S. Benk, M. Tarnawski, J. P. Spatz and R. Dimova, ACS Nano, 2023, 17, 11957–11968.
50. R. Jahn and R. H. Scheller, Nat. Rev. Mol. Cell Biol., 2006, 7, 631–643.
51. S. Koike and R. Jahn, Biochem. J., 2022, 479, 273–288.
52. D. Mion, L. Bunel, P. Heo and F. Pincet, FEBS Open Bio, 2022, 12, 1958–1979.
53. R. Lipowsky, R. Ghosh, V. Satarifard, A. Sreekumari, M. Zamaletdinov, B. Różycki, M. Miettinen and A. Grafmüller, Biomolecules, 2023, 13, 926.
54. K. Wang, M. E. Villanueva, F. Caporaletti, R. P. White, J. E. G. Lipson, S. Napolitano and P. Losada-Pérez, Small, 2026, 22, e12844.
55. L. Martnez, R. Andrade, E. G. Birgin and J. M. Martnez, J. Comput. Chem., 2009, 30, 2157–2164.
56. A. I. Jewett, D. Stelter, J. Lambert, S. M. Saladi, O. M. Roscioni, M. Ricci, L. Autin, M. Maritan, S. M. Bashusqeh and T. T. Keyes, et al., J. Mol. Biol., 2021, 433, 166841.
57. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. Crozier, P. J. I. Veld, A. Kohlmeyer, S. G. Moore and T. D. Nguyen, et al., Comput. Phys. Commun., 2022, 271, 108171.
58. M. M. Anila, R. Ghosh and B. Różycki, Soft Matter, 2023, 19, 3723–3732.

---