The SNARE proteins syntaxin-1, SNAP-25 and synaptobrevin mediate neurotransmitter release by forming tight SNARE complexes that fuse synaptic vesicles with the plasma membranes in microseconds. Membrane fusion is generally explained by the action of proteins on macroscopic membrane properties such as curvature, elastic modulus and tension, and a widespread model envisions that the SNARE motifs, juxtamembrane linkers and C-terminal transmembrane regions of synaptobrevin and syntaxin-1 form continuous helices that act mechanically as semi-rigid rods, squeezing the membranes together as they assemble (‘zipper’) from the N- to the C-termini. However, the mechanism underlying fast SNARE-induced membrane fusion remains unknown. We have used all-atom molecular dynamics simulations to investigate this mechanism. Our results need to be interpreted with caution because of the limited number and length of the simulations, but they suggest a model of membrane fusion that has a natural physicochemical basis, emphasizes local molecular events over general membrane properties, and explains extensive experimental data. In this model, the central event that initiates fast (microsecond scale) membrane fusion occurs when the SNARE helices zipper into the juxtamembrane linkers which, together with the adjacent transmembrane regions, promote encounters of acyl chains from both bilayers at the polar interface. The resulting hydrophobic nucleus rapidly expands into stalk-like structures that gradually progress to form a fusion pore, aided by the SNARE transmembrane regions and without clearly discernible intermediates. The propensity of polyunsaturated lipids to participate in encounters that initiate fusion suggests that these lipids may be important for the high speed of neurotransmitter release.

The methodology used in system setups, equilibration and production MD simulations was analogous to that described previously (1). The vesicle generated previously (1) was adapted for each system by moving lipids manually to accommodate different positions of the SNARE TM regions or their absence. A flat square bilayer of 30.5 nm x 30.5 nm was built at the CHARMM-GUI website (2) (https://charmm-gui.org/) for the fusion2g system and was adapted for all other systems (also moving lipids manually to accommodate SNARE TM regions) except for the snfreet350 system, which used a smaller bilayer adapted from the qscv simulation of ref. (1). The lipid compositions of the original vesicle and flat bilayers are described in Table 1 of the manuscript, and those of each system were only slightly altered by adding or removing a few lipids as needed. The number of atoms, box dimensions and simulation temperatures of each system are also listed in Table 1 of the manuscript. All systems were solvated with explicit water molecules (TIP3P model), adding potassium and chloride ions as needed to reach a concentration of 145 mM and make the system neutral.

               All systems were energy minimized using double precision, heated to the desired temperature over the course of a 1 ns MD simulation in the NVT ensemble with 1 fs steps, and equilibrated to 1 atm for 1 ns in the NPT ensemble using isotropic Parrinello-Rahman pressure coupling and 2 fs steps (3). NPT production level MD simulations were performed for the times indicated in Table 1 using 2 fs steps, isotropic Parrinello-Rahman pressure coupling and a 1.1 nm cutoff for non-bonding interactions. Nose-Hoover temperature coupling (4) was used separately for three groups: i) protein atoms; ii) lipid atoms; and iii) water and KCL. Periodic boundary conditions were imposed with Particle Mesh Ewald (PME) (5) summation for long-range electrostatics. The default mixed precision was used in all MD simulations. The speeds of the production simulations ran on Frontera at TACC were typically about 28 ns/day for a typical system of about 5.2 million atoms using 40 nodes.

(1) J. Rizo, L. Sari, Y. Qi, W. Im, M. M. Lin, All-atom molecular dynamics simulations of Synaptotagmin-SNARE-complexin complexes bridging a vesicle and a flat lipid bilayer. Elife 11, e76356 (2022).

(2) S. Jo, T. Kim, V. G. Iyer, W. Im, CHARMM-GUI: a web-based graphical user interface for CHARMM. J Comput Chem 29, 1859-1865 (2008).

(3) Parrinello M, Rahman A. 1981. Polymorphic transitions in single crystals: A new molecular dynamics method. Journal of Applied Physics 52:7182–7190

(4) Hoover T. 1985. Canonical dynamics: Equilibrium phase-space distributions. Physical Review. A, General Physics 31:1695–1697.

(5) Darden T, York D, Pedersen L. 1993. Particle mesh Ewald: an N·log(N) method for Ewald sums in large systems. The Journal of Chemical Physics 98:10089–10092.