Nanoparticle-induced biomembrane fusion: unraveling the effect of core size on stalk formation

Abstract

Membrane fusion in vitro is a strategy to load model or cell-derived vesicles with proteins, drugs, and genetic materials for theranostic applications. It is thus crucial to develop strategies to control the fusion process, also through synthetic fusogenic agents. Ligand-protected, membrane-penetrating gold nanoparticles (Au NPs) can facilitate membrane fusion, but the molecular mechanisms remain unresolved. Here, we tackle NP-induced stalk formation using a coarse-grained molecular dynamics approach and enhanced sampling techniques. We show that smaller (2 nm in diameter) NPs lead to a lower free energy barrier and higher stalk stability than larger NPs (4 nm). We demonstrate that this difference is due to a different ligand conformational freedom, which in turn depends on the Au core curvature. Our study provides precious insights into the mechanisms underlying NP-mediated membrane fusion, while our computational approach is general and applicable to studying stalk formation caused by other fusogenic agents.

---

Membrane fusion is fundamental for many biological activities such as membrane trafficking,¹ synaptic transmission,^2,3 fertilization,^4,5 and viral infection.⁶ It is a complex process involving the concerted movement and rearrangement of many different molecules (mainly lipids and proteins). The most likely fusion path, from two approaching bilayers to the opening of the fusion pore, includes intermediate metastable minima (Fig. 1a). At first, the two separated membranes must be brought close to each other, typically at <5 nm distance. Then, a point-like protrusion of a lipid tail can initiate the hydrophobic contact between the membranes, leading to the so-called stalk, a configuration allowing the mixing of the lipids from the outer leaflets.^7–10 Stalk elongation may lead to a second metastable configuration, the hemifusion diaphragm, and then the process is completed with the opening and expansion of a fusion pore.^10–12 The fusion pore opening could directly evolve from the stalk configuration.¹² The different minima characterizing this complex multi-step process can be separated by relatively high free energy barriers.^8,9,13

Fig. 1 Fusion intermediates and definition of the chain coordinate. (a) Scheme of the typical free energy landscape along the fusion path. It exhibits three free energy barriers: one for the stalk formation, one for the expansion of the hemifusion diaphragm, and one for the opening of the fusion pore. (b) A scheme of the cylinder used by Hub and co-workers to define the chain coordinate ξ_ch.⁹ The cylinder is divided into slices and placed in between two adjacent membranes. When the membranes are distinct, most slices are empty of lipid tails, while in the stalk configuration, all slices are filled.

In vivo, overcoming the fusion barriers is generally achieved thanks to the presence of specialized proteins, such as the SNARE complex, which catalyze the fusion of synaptic vesicles with the plasma membrane.^12,14In vitro, achieving controlled fusion by means of artificial fusogens would be critical to realize artificial vesicles with controlled membrane composition and cargo.^15,16 Gold nanoparticles (Au NPs), functionalized with anionic (11-mercapto-1-undecanesulfonate – MUS) and neutral, hydrophobic (octanethiol – OT) ligands, are promising synthetic fusogenic agents. Indeed, a recent study by Tahir et al. demonstrated that MUS:OT NPs induce fusion between pure phosphatidylcholines (PC) liposomes,¹⁷ and our group showed that the liposome cholesterol content can finely tune MUS:OT NP-induced fusion events.¹⁸

Molecular Dynamics (MD) simulations, thanks to their high spatiotemporal resolution, have been successfully employed to shed light on the molecular mechanisms involved in the different stages of the fusion process, either spontaneous^19–22 or mediated by fusion agents.^17,23–25 The presence of high free energy barriers makes the use of atomistic models challenging and calls for coarse-grained models and/or enhanced sampling techniques. Thanks to coarse-grained modelling, spontaneous NP-induced stalk formation has been observed in unbiased MD, but only increasing the temperature much above the physiological one.¹⁸ Therefore, the derivation of the free energy profile of stalk formation requires enhanced sampling, and the challenge is the definition of a collective variable (CV) that approximates well the reaction coordinate. Finding a CV suitable to describe stalk formation is arduous because of the high number of molecules involved and their high conformational freedom, and a poor CV could easily lead to hysteresis in the free energy calculation along the stalk-formation and stalk-destruction pathways.²⁶

The group of J. S. Hub recently solved these issues by proposing the so-called “chain coordinate” ξ_ch,^9,26–28 designed to study stalk formation between two flat lipid bilayers and implemented directly into the Gromacs MD suite. In order to quantify the degree of connectivity between the hydrophobic cores of two lipid membranes separated by a water layer, one can imagine to track the filling/unfilling of a cylindrical volume placed in between the two membranes, with its axis aligned to the membrane normal, as shown in Fig. 1b. ξ_ch represents the cylinder's degree of filling by hydrophobic lipid tails. In the stalk configuration, all slices are filled with hydrophobic lipid tails, and ξ_ch approaches unity, while when the membranes are separated, several slices are empty, and ξ_ch is ≪ 1 (see ESI† for further details on the ξ_ch original definition). Hub and co-workers exploited ξ_ch in umbrella sampling calculations and showed how lipid membrane composition impacts the stalk formation free energy path.⁹

Here we build on the original definition of the ξ_ch chain coordinate, modifying it to account for the presence of the fusogenic agent, a MUS:OT Au NP embedded in one of the two lipid bilayers. Our previous MD simulations based on a coarse-grained force field revealed that the stalk rim constantly forms over a single Au NP,¹⁸ consequently to the formation of a hydrophobic contact between a lipid tail of the top bilayer and transient hydrophobic defects on the NP, water-facing surface. This mechanism is consistent with the one described in the literature for NP-free membranes, with the only difference being that the contact is a lipid-NP one and not a lipid–lipid one.⁹ Therefore, we decided to redefine the chain coordinate ξ_ch to make it suitable to deal with a NP-induced stalk formation (Fig. 2), implementing the following two changes.

Fig. 2 A new implementation of the chain coordinate. The cylinder is now defined with respect to the NP, following it in the xy plane and being located along z at a distance h from the NP center of mass. (a) When the membranes are still distinct, the cylinder is partially empty (ξ_ch ≪ 1). (b) When a stalk is formed over the NP, the cylinder is almost full of hydrophobic beads (ξ_ch ∼ 1). Lipid heads are in orange, lipid tails in cyan, the NP core in yellow, MUS ligands in red, OT ligands and cholesterol are not shown for clarity.

First, to have the cylinder always follow the NP diffusing in the membrane plane, we set the in-plane x and y coordinates of the cylinder axis to coincide with the x and y of the NP center of mass (COM). Second, to account for the system asymmetry along the membrane normal and the possible differences in NP size, we set the z of the cylinder center of mass to a user-defined distance h from the NP COM z coordinate (see ESI† for more details). As shown in Fig. 2, even in the presence of NP diffusion and oscillation along the z direction, the cylinder is automatically centered over the NP, thus monitoring the degree of hydrophobic connection established between the two lipid bilayers just on top of the NP.

To test our newly defined CV, we set up a system consisting of two DOPC bilayers with 30% mol cholesterol, a lateral size of ∼19 nm, and a single NP with a 2 nm core fully embedded in the lower bilayer (see ESI Fig. S1a†). After minimization of our solvated system, we equilibrated it and then derived stalk formation free energy profiles along ξ_ch using the Umbrella Sampling (US) technique. We extracted the starting configuration of each of the 16 windows from a simulation in which we slowly pulled the system along ξ_ch to drive stalk formation (see ESI† for more details). We set the shift parameter h to 2.75 nm to account for the size of the NP core. The free energy profile (reported in Fig. 3a) was constructed from the umbrella histograms with the weighted histogram analysis method (WHAM),²⁹ as implemented in the gmx wham Gromacs tool, and the error was calculated from a bootstrap procedure.³⁰ The free energy profile, shown in red in Fig. 3a, shows a local minimum at ξ_ch = 0.31 that corresponds to the state in which the two membranes are separated (pre-stalk state). The second minimum, at ξ_ch = 0.94, corresponds to the state in which the stalk is fully formed. In this case, the stalk is thermodynamically favored, being lower in free energy of about 30 kJ mol⁻¹, and the barrier to stalk formation is about 20 kJ mol⁻¹. This picture is consistent with previous unbiased simulations, which were run at a higher temperature and estimated a barrier of about 15 kJ mol⁻¹ in the presence of cholesterol.¹⁸ Our calculation here demonstrates that, at physiological temperature, the stalk formed on top of the NP is highly stable, with a barrier to stalk destruction larger than 50 kJ mol⁻¹.

Fig. 3 Free energy profiles of stalk formation. (a) Free energy profiles of stalk formation (red) and destruction (orange) as a function of ξ_ch for NP with a core size of 2 nm. The absence of hysteresis demonstrates the reliability of the CV. (b) Free energy profiles of stalk formation for a 2 nm core (red) and a 4 nm core (blue). The barrier to stalk formation in the second case is higher, while the stalk becomes much less stable, resulting in a minimum with the same free energy as the apposed membranes one.

To test the consistency of our method and exclude the presence of significant hysteresis, we calculated the same free energy profile using initial configurations taken from a reverse pulling simulation, in which the system has been brought from the stalk state back to the pre-stalk state. We refer to this second profile as the “profile of stalk destruction” (Fig. 3a, orange). In the ideal case, if the free energy profiles are at convergence, the free energy profiles of stalk formation and destruction should be identical. Comparing the two profiles in Fig. 3a, it can be noticed that they overlap within the error bars, proving the absence of hysteresis and thus the reliability of our CV.

Once verified the appropriateness of our method, we tackled the effect of the NP core size on the stalk formation process. To this purpose, we constructed a second system, still containing two membranes composed of 70% DOPC and 30% cholesterol, but with a single, embedded NP with a gold core diameter of 4 nm (see ESI Fig. S1b†). Setting up this new system, we set the lipid hydration to the same value as in the previous simulations. The free energy profile of stalk formation was calculated as in the previous case (see ESI† for more details), with the only different parameter being the h shift, which was increased to 3.75 nm to account for the new core radius. Once again, we verified the absence of significant hysteresis by calculating the free energy profile of stalk destruction, and obtained an almost perfect overlap with the one of stalk formation (see ESI Fig. S2†).

Fig. 3b compares the free energy profiles of stalk formation in the two systems. We can notice a different barrier and a different depth of the stalk minimum. In the presence of the bigger NP, the stalk is not thermodynamically favored anymore, becoming equally probable to the pre-stalk state. Furthermore, the barrier to stalk formation is higher, suggesting a slower kinetics of stalk formation for the large NPs. This result is consistent with recent experimental observations in which MUS:OT NPs with a core size of ∼2.5 nm promoted hemifusion between adjacent bilayers, while larger NPs exhibited only surface binding to membranes.³¹ To further prove the size effect, we built a CG model of an intermediate-size nanoparticle (3 nm of core diameter) and calculated the associated free energy profile of stalk formation, following the same protocol adopted for the other two cases. The result, reported in ESI Fig. S6,† confirms the size effect: both the free energy barrier and the depth of the stalk minimum are intermediate between those of the 2 nm case and those of the 4 nm case. All the results obtained at three different sizes are now resumed in Table 1.

Table 1 A summary of the results obtained on the barrier to stalk formation and the depth of the stalk free energy minimum (which are both calculated with respect to the pre-stalk state)

NP core diameter (nm)	Height of the barrier (kJ mol^−1)	Stalk minimum depth (kJ mol^−1)
2	19.5 ± 1	−33 ± 2
3	22 ± 1	−8.5 ± 1
4	26 ± 1	0 ± 1

---

In order to analyze and explain the molecular origin of this size effect, we calculated and compared a few crucial quantities in the two extreme cases (2 and 4 nm), relating them to the differences observed in the free energy profiles of Fig. 3b. The molecular trigger for stalk formation is the establishing of a first contact between the hydrophobic tail of a lipid protruding from the facing bilayer (in our set-up, the one above the NP) and a transient hydrophobic defect in the NP ligand shell that is exposed to the water phase (see ESI Fig. S3†). This is a rare event, whose rate depends on the height of the free energy barrier separating the pre-stalk state from the stalk state.

We wanted to test the hypothesis that the more hydrophobic surface is exposed by the NP to the aqueous environment, the more likely it is that a hydrophobic NP-lipid contact is established. To quantify this, we defined the H-LASA as the hydrophobic lipid-accessible-surface-area of the NP divided by the total lipid accessible surface area of the NP (see ESI† for more details). We then monitored the time evolution of the H-LASA during a simulation (see ESI Fig. S4†) and recorded its maximum value.

Fig. 4a shows the maximum H-LASA obtained for NPs of both sizes when embedded in the membrane in the pre-stalk configuration. The maximum H-LASA anti-correlates with the barrier to stalk formation shown in Fig. 3, being significantly larger for the 2 nm core size NP than for the 4 nm one. This difference can be explained with simple geometric considerations. As the core radius becomes bigger, the average surface curvature decreases, while the ligand grafting density and ligand length remain the same. Therefore, the interface of the large NP will have both a higher ligand density and a higher surface charge density than a smaller NP. Both the ligand density and the charge density limit the ligand conformational freedom and make the formation of hydrophobic patches more difficult. The same effect is present when the NPs are solvated in water (see ESI Fig. S5†), and it is responsible for hindering the hydrophobic aggregation between the large NPs.³²

Fig. 4 Microscopic interpretation of the NP core size effect on stalk formation. (a) Maximum of H-LASA in the pre-stalk configuration for the NPs with a diameter of 2 nm and 4 nm. (b) Number of hydrophobic contacts in the stalk state for the two NPs. (c) Number of hydrophobic lipid beads present in the cylinder in the stalk state for the two different size NPs. All the error bars in panels (a), (b), and (c) are standard errors calculated from simple block averages. (d) Representative snapshots with views from the side and from the top of the NPs of the two sizes in the stalk state. The opening of the MUS-charged terminals due to the formation of the stalk above the nanoparticle is highlighted. We show only the MUS charged terminals (in transparent red) and the NPs' core (in yellow). The width of the “hole” corresponds to the diameter of the stalk above the NP, which is different in the two cases.

The impact of the NP core size on the stalk stability is even more pronounced than the effect on the barrier to stalk formation. What is the molecular origin of the higher stalk stability in the case of the small NP? A stalk is a hydrophobic connection between two membranes, and we expect its stability to be correlated with the stalk geometry and composition. A larger/denser stalk is expected to provide a more stable hydrophobic connection between the membranes, even when mediated by the NP. Are the stalk dimensions and densities different in the two systems investigated here?

We calculated the number of hydrophobic contacts between the NP and the lipids in the stalk region. More precisely, we calculated the variation of the number of contacts between the hydrophobic beads of the NP ligands and those of the membranes between the pre-stalk and the stalk state. The results, shown in Fig. 4b, correlate well with the different stalk stability. In the case of the 2 nm NP, the system takes advantage of the possibility of forming more NP-lipid hydrophobic contacts and thus manages to create a more stable stalk.

To get a quantitative estimation of the stalk diameters, we counted all the hydrophobic beads contained in the 18 slices of the cylinder closer to the NP and thus located in the central part of the stalk region (see ESI† for more details). The results are shown in Fig. 4c. The stalk contains more hydrophobic beads when the NP is smaller. From the values obtained, being the cylinder height constant, it is straightforward to derive the ratio between the average diameter of the stalk in the case of a 4 nm core and the one in the case of a 2 nm core, which turns out to be about 0.7. The difference in stalk diameter between the two systems becomes evident in Fig. 4d, where we show representative snapshots of the NPs in the stalk state. The opening of the charged MUS terminals due to the formations of the stalk above the NP is highlighted by our representation, where the diameter of the “hole” corresponds to the diameter of the hydrophobic connection between the NP and the membrane above. The hole (and thus the stalk) is larger in the system containing the small NP. Again, the molecular origin of this difference is the higher conformational freedom of MUS ligands on the small NP.

Conclusions

In conclusion, we developed a reliable procedure to calculate the free energy of stalk formation induced by functionalized NPs. Our approach was based on adapting the chain coordinate⁹ to the presence of membrane-embedded fusogenic agents. We calculated the free energy profiles of stalk formation via umbrella sampling, revealing that the NP core size can significantly impact both the kinetics and the thermodynamics of the process. A smaller NP core leads to higher conformational freedom of charged ligands, which can more easily rearrange to allow for the first hydrophobic contact between the NP and the facing lipids, a contact that triggers stalk formation. Coherently, when the core is smaller, the open-ligand configuration is more stable, and it allows for the formation of a larger, durable stalk, driving its thermodynamic stability over the pre-stalk state.

Stalk formation is a fundamental step of membrane fusion, and thus our results have important implications for the design of artificial fusogenic agents. In particular, we highlight that ligand flexibility and conformational freedom are critical features that can be exploited to increase the efficiency of the fusion process. The biased protocol developed here will be used to study the influence of the membrane composition on NP-induced stalk formation and will be combined with further in silico studies to reveal how the system can go beyond the stalk state, completing the fusion process. It is worth noting that our collective variable definition is quite general and, thus, not limited to the particular type of NP studied herein. Indeed, our protocol could be easily adapted to stalk formation induced by different artificial or biological agents, helping develop new strategies to induce and control fusion events.

Author contributions

Giorgia Brosio: investigation, methodology, software, writing – first draft. Giulia Rossi: methodology, supervision, resources, writing – review & editing. Davide Bochicchio: conceptualization, methodology, supervision, writing – first draft, writing – review & editing.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors acknowledge funding by MIUR – DIFI Dipartimento di Eccellenza 2018–2022 for computational resources.

Notes and references

1. R. Jahn and T. C. Südhof, Annu. Rev. Biochem., 1999, 68, 863–911.
2. Z. Yang, L. Gou, S. Chen, N. Li, S. Zhang and L. Zhang, Front. Mol. Neurosci., 2017, 10, 168.
3. T. Söllner and J. E. Rothman, Trends Neurosci., 1994, 17, 344–348.
4. P. Primakoff and D. G. Myles, FEBS Lett., 2007, 581, 2174–2180.
5. B. M. Gadella and J. P. Evans, Adv. Exp. Med. Biol., 2011, 713, 65–80.
6. S. C. Harrison, Nat. Struct. Mol. Biol., 2008, 15, 690–698.
7. M. M. Kozlov, S. L. Leikin, L. V. Chernomordik, V. S. Markin and Y. A. Chizmadzhev, Eur. Biophys. J., 1989, 17, 121–129.
8. Y. Kozlovsky and M. M. Kozlov, Biophys. J., 2002, 82, 882–895.
9. C. S. Poojari, K. C. Scherer and J. S. Hub, Nat. Commun., 2021, 12, 6594.
10. L. V Chernomordik and M. M. Kozlov, Nat. Struct. Mol. Biol., 2008, 15, 675–683.
11. H. J. Risselada, G. Bubnis and H. Grubmüller, Proc. Natl. Acad. Sci., 2014, 111, 11043–11048.
12. T. Lu and H. Guo, Adv. Theory Simul., 2019, 2, 1900032.
13. A. Grafmüller, J. Shillcock and R. Lipowsky, Biophys. J., 2009, 96, 2658–2675.
14. R. Jahn and R. H. Scheller, Nat. Rev. Mol. Cell Biol., 2006, 7, 631–643.
15. M. Piffoux, A. K. A. Silva, C. Wilhelm, F. Gazeau and D. Tareste, ACS Nano, 2018, 12, 6830–6842.
16. D. Mukherjee, D. Paul, S. Sarker, Md. N. Hasan, R. Ghosh, S. E. Prasad, P. K. Vemula, R. Das, A. Adhikary, S. K. Pal and T. Rakshit, ACS Appl. Bio Mater., 2021, 4, 8259–8266.
17. M. A. Tahir, Z. P. Guven, L. R. Arriaga, B. Tinao, Y.-S. S. Yang, A. Bekdemir, J. T. Martin, A. N. Bhanji, D. Irvine, F. Stellacci and A. Alexander-Katz, Proc. Natl. Acad. Sci., 2020, 117, 18470–18476.
18. E. Canepa, D. Bochicchio, G. Brosio, P. H. J. Silva, F. Stellacci, S. Dante, G. Rossi and A. Relini, Small, 2023, 2207125.
19. H. J. Risselada, C. Kutzner and H. Grubmüller, ChemBioChem, 2011, 12, 1049–1055.
20. Y. G. Smirnova, S.-J. Marrink, R. Lipowsky and V. Knecht, J. Am. Chem. Soc., 2010, 132, 6710–6718.
21. S. Kawamoto and W. Shinoda, Soft Matter, 2014, 10, 3048–3054.
22. M. Schick, J. Stat. Phys., 2011, 142, 1317–1323.
23. S. Baoukina and D. P. Tieleman, Biophys. J., 2010, 99, 2134–2142.
24. R. M. Bhaskara, S. M. Linker, M. Vögele, J. Köfinger and G. Hummer, ACS Nano, 2017, 11, 1273–1280.
25. R. Santiago and R. Reigada, Biochim. Biophys. Acta, Gen. Subj., 2019, 1863, 723–731.
26. N. Awasthi and J. S. Hub, J. Chem. Theory Comput., 2016, 12, 3261–3269.
27. W. F. D. Bennett and D. P. Tieleman, J. Chem. Theory Comput., 2011, 7, 2981–2988.
28. C. L. Ting, N. Awasthi, M. Müller and J. S. Hub, Phys. Rev. Lett., 2018, 120, 128103.
29. S. Kumar, J. M. Rosenberg, D. Bouzida, R. H. Swendsen and P. A. Kollman, J. Comput. Chem., 1992, 13, 1011–1021.
30. B. Efron, Ann. Math. Stat., 1979, 7, 1–26.
31. P. U. Atukorale, Z. P. Guven, A. Bekdemir, R. P. Carney, R. C. Van Lehn, D. S. Yun, P. H. Jacob Silva, D. Demurtas, Y.-S. Yang, A. Alexander-Katz, F. Stellacci and D. J. Irvine, Bioconjugate Chem., 2018, 29, 1131–1140.
32. E. Lavagna, D. Bochicchio, A. L. De Marco, Z. P. Güven, F. Stellacci and G. Rossi, Nanoscale, 2022, 14, 6912–6921.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3na00430a

---