## Abstract

Neurotransmitter release from neuronal terminals is governed by synaptic vesicle fusion. Vesicles filled with transmitters are docked at the neuronal membrane by means of the SNARE machinery. After a series of events leading up to the fusion pore formation, neurotransmitters are released into the synaptic cleft. In this paper, we study the mechanics of the docking process. A continuum model is used to determine the deformation of a spherical vesicle and a plasma membrane, under the influence of SNARE-machinery forces and electrostatic repulsion. Our analysis provides information on the variation of in-plane stress in the membranes, which is known to affect fusion. Also, a simple model is proposed to study hemifusion.

## 1. Introduction

Neuronal transmitters are packed into vesicles and released from synaptic terminals by fusion of the vesicles with the plasma membrane. Synaptic vesicle fusion is mediated by the protein complex termed the SNARE. Details of the mechanism are still unclear, because it occurs at the microsecond scale, and thus it is difficult to capture the process using current imaging techniques. The release of transmitters into the synaptic cleft involves two major steps: (i) docking of vesicle on the plasma membrane and (ii) membrane fusion and subsequently pore opening. The central role in this entire mechanism is played by a specialized group of proteins termed the SNARE (soluble N-ethylmaleimide-sensitive-factor attachment protein receptor) complex [1–4]. More specifically, the process involves two transmembrane proteins, synaptobrevin 2 (Syb) and syntaxin 1 (Syx), and the protein SNAP-25 [5]. Syb is attached at one end to the vesicle by embedding a hydrophobic transmembrane domain [6] into the vesicle membrane. It is hence referred to as ‘v-SNARE’. The protein Syx embeds a transmembrane domain [7], and the protein SNAP-25 is anchored with a palmitoyl chain [8] onto the plasma membrane. The proteins attached to the target membrane are termed as ‘t-SNAREs'. When a vesicle is near the plasma membrane, t-SNAREs form an acceptor site for the v-SNARE leading to the formation of the SNARE complex [9,10]. This SNARE complex consists of four helices (contributed by the Syb, Syx and SNAP-25), zippered into a tight bundle. It has been suggested that several SNARE complexes are involved in the docking process and the zippering action of these SNAREs provides a force to counter the repulsive electrostatic force between the membranes. Gao *et al.* [11] using optical tweezers controlled pulling experiments, found that the SNARE-machinery generates forces of an order of 2–20 pN, depending on the level of zippering. Recently, calculations of the electrostatic repulsion between the vesicle and the membrane suggested that the repulsive force changes strongly with distance between the vesicle and membrane and is of the order of 100–200 pN for a separation of 1 nm between a typical synaptic vesicle and plasma membrane [12].

Synaptic vesicle and plasma membranes are bilayers made up of relatively long amphipathic lipid molecules, with their hydrophilic heads in the aqueous solution on either side of their bilayer and the hydrophobic tails in the interior. Lipid membranes are generally regarded as two-dimensional fluids that, from a modelling point of view, conserve their area during deformation [13,14]. Because the outer surface of the bilayer usually carries a net negative charge, considerable force is required to overcome the electrostatic repulsion between the vesicle and plasma membranes to dock the former. It has been hypothesized that the energy released during the SNARE zippering is used to overcome the energy barrier of the electrostatic interaction [15].

Coarse-grained continuum models for the contact mechanics of vesicle membranes have been studied by various authors [16–28]. Recently, Blount *et al.* [29] analysed the problem of a pressurized cylindrical vesicle interacting with a rigid substrate under a potential which has short-range repulsion and long-range attraction. Our approach in this work is similar, except that our vesicle is spherical, our substrate is a lipid membrane and the attractive potential is replaced by the zipping force exerted by the SNARE-machinery, which is modelled as a concentrated line load acting on both membranes. Electrostatic repulsion between the membranes is modelled using the Debye–Hückel (DH) theory [12].

## 2. Geometry and model

The geometry is shown schematically in figure 1. We use a cylindrical coordinate system (*r*, *θ*, *z*) with *θ* the angle of revolution about *z*-axis. To simplify the calculation, we assume axisymmetry; that is, we model the docking process by prescribing a circle of line forces of magnitude *F* on a spherical vesicle of radius *R* (figure 1) as well as on the plasma membrane. These forces represent the zipping of the SNARE-machinery and counter the repulsive electrostatic forces between them. As shown in figure 1, the line force acts along a latitude of the undeformed vesicle and is constrained to remain normal to the deformed surface. The location of the latitude is specified by the arc length *S*_{0} of a cross section in the reference configuration, which is taken to be a spherical vesicle. Because the plasma membrane is very large compared with the vesicle radius, its reference configuration is taken to be a flat circular membrane of radius *L* under pretension *T*_{0}. We allow two different types of line forces acting on the plasma membrane. Both sets of forces act on a circle of radius *S*_{0} and have the same magnitude *F*. The first set is assumed to be always normal to the deformed plasma membrane, whereas the second set is always directed opposite to the force on the vesicle (figure 1).

The double layer model of charged surfaces in electrolytes has been successfully used to model the electrostatics near the lipid bilayers [30,31]. In this work, we follow the double layer model used by Bykhovskaia *et al.* [12] where electrostatic repulsion between the vesicle and plasma membrane is determined by the solution of the DH equation. McIntosh *et al.* [32] observed that when the gap between the membranes is greater than 1 nm, the DH equation predicts membrane interactions consistent with their experiments. For membrane separation below 1 nm, which is close to the Debye length *l*_{D} ≈ (0.67 nm) in our system, hydration repulsion, van der Waals attraction and nonlinearity in the electrostatics come into play. Thus, as long as the gap between the membranes is greater than the Debye length, the solution of the DH equation should give a reasonable description of the electrostatics. It will be shown in the Results section that, for the issues addressed in this paper, the separation between the membranes is always greater than twice the Debye length, ensuring that our electrostatics is consistent with the DH approximation. In our model, the electrostatic interaction between curve surfaces is calculated by the Derjaguin approximation [33], which assumes that locally the surfaces are flat, so interaction between two material points on the different membranes separated by a distance of *δ* (figure 1*b*) can be described by determining the force per unit area between two infinite parallel planes separated by the same distance. The DH equation, appropriate for this geometry is given by equation (2.1*a*). The repulsive force per unit area between the two planes, *F*_{e}, is obtained by solving equation (2.1*a*) using the constant charge boundary conditions (equation (2.1*b*)). Details of solution are given in the electronic supplementary material, eq. 105.2.1a–cwhere *σ*_{ν} (=−0.025 C m^{−2}) and *σ*_{pm} (=−0.070 C m^{−2}) are the surface charge densities of the vesicle and plasma membrane, respectively, *l*_{D} (=0.67 nm) is the Debye length, *ɛ* (=80) is the relative permittivity of water and *ɛ*_{0} is the permittivity of vacuum. Our choice of surface charge densities above is based on the work of Bykhovskaia *et al.* [12]. It is about two to three times higher than those reported by Pekker & Shneider [34]. However, as noted later (figure 6), the tension profile between membranes is insensitive to the magnitude of charge densities within this range. We also check that the equilibrium distances between the membranes are also not very sensitive to the magnitude of charge densities within this range (see the electronic supplementary material, figure S5).

The membranes are deformed by the zipping force of the SNARE complexes and electrostatic repulsion force. We model this deformation using continuum theory. The strain energy densities *W* of both membranes are given by
2.2where *H* is the mean curvature and *c* ≈ 10–20*k*_{B}*T* is the bending rigidity of the lipid bilayer [35].

### 2.1. Governing equations for vesicle membrane

In the following, we use the formulation of Jenkins [17,36] and Long *et al.* [37] to derive the governing equations for the deformation of the spherical vesicle and the flat plasma membrane. The undeformed configuration of the vesicle is a sphere of radius *R* with arc length in a cross section denoted by *S*, whereas the plasma membrane occupies the interior of a circle of radius . We introduce the notation *ϕ* to denote the angle made by the tangent to a point on the cross section of the deformed membrane in the (*r*, *z*) plane with the *z*-axis (figure 2*a*). The osmotic pressure of the vesicle is denoted by *p*_{0}. The forces acting on the vesicle and plasma membrane due to SNARE complexes (several SNARE complexes can be attached to a vesicle, and the forces due to these vesicles are distributed uniformly on a closed circular arc on the vesicle) are denoted by *F* and –*F*, respectively (figure 1). Forces are resolved into a normal (*F*_{n}) and tangent component (*F*_{t}) with respect to the deformed vesicle and plasma membrane. The electrostatic interaction force *F*_{e} is given by equation (2.1) and is assumed to act in the *z*-direction.

The equations describing the deformation involve the angle *ϕ*, the mean curvature *H*, the deformed arc length *ξ*, the deformed coordinates of a generic material point (*r*, *z*) which has an arc-length coordinate *S* in the undeformed configuration. The force variables relevant to the calculation involve the shear force *Q* and *d* an integration constant which determines the tension *T* in the membranes.

To reduce the number of parameters in our simulations, we introduce normalized variables which are indicated by a horizontal bar. All distances are normalized by *R*, the radius of the vesicle. As the bending rigidity *c* has units of energy, we use it to normalize force per unit length quantities, i.e. the out of plane shear *Q* and in-plane tension *T* are normalized by *c*/*R*^{2}. Also, force per unit area quantities, *p*_{0}, *F*_{e}, *F*_{t} and *F*_{n} are made dimensionless by dividing with *c*/*R*^{3}. These variables are summarized as follows:2.3In all simulations, we set the tangential component of SNARE force *F*_{t} on the vesicle to zero. In the electronic supplementary material, we derived in detail the six ordinary differential equations governing the deformation of the vesicle membrane (the final forms of these equations are given in the electronic supplementary material, eqs. 85 and 87):2.4a–fwhere the dot denotes differentiation with respect to the normalized undeformed arc-length , and2.4gEquation (2.4*a*) represents force balance in the normal direction (electronic supplementary material, eq. 83). Equations (2.4*b*,*c*) are rearranged forms of the mean curvature (electronic supplementary material, eq. 74) and shear force (electronic supplementary material, eq. 77) definitions, respectively. Equations (2.4*d*,*e*) represent the geometric relationship between the variables (electronic supplementary material, eq. 65), whereas equation (2.4*f*) is obtained from the force balance in tangential direction (electronic supplementary material, eq. 84). The generalized pressure in equation (2.4*a*) is related to the normalized osmotic pressure , the electrostatic force per unit area and the normal component of the line load applied at , by2.4hwhere is the Dirac delta function.

These differential equations are supplemented with the boundary conditions2.5a–f

The boundary conditions defined above describe the symmetry in the vesicle geometry. About the symmetry (*z*) axis, the curve has zero slope and the out of plane shear is zero, at both and *π*. Also, for the continuity of the geometry, we impose at both and *π*.

The notation for positive shear force and tension is described in figure 2*b*. Finally, the expression for the normalized in-plane tension *T* in both the vesicle and plasma membranes is given by2.6

### 2.2. Governing equations for plasma membrane

The governing equations for the deformation of the plasma membrane are the same as equations (2.4*a–f*) (see the electronic supplementary material for details) except that equation (2.4*g*) must be replaced by2.7

This change is due to the difference between the reference configurations (one is a sphere and the other a flat surface). The boundary conditions are2.8a–f

The first three boundary conditions equations (2.8*a–c*) are due to axisymmetry. Equation (2.8*f*) states that the tension in the plasma membrane approaches the pretension at the boundary. This boundary condition along with equations (2.8*d*,*f*) allow the plasma membrane to deflect only in horizontal direction. Had we replaced this boundary condition with a clamped condition, the deflection everywhere would be zero because of area incompressibility.

The coupled differential equations given by equations (2.4*a*–*g*) and equation (2.7) with the boundary conditions given by equations (2.5*a–f*) and (2.8*a–f*) are solved using the Matlab bvp4c solver. The input parameters for the solver are the osmotic pressure *p*_{0} across the vesicle membrane, which remains fixed throughout the deformation, SNARE-machinery force parameters (*S*_{0} and *F*), electrostatic force (*F*_{e}) and pretension (*T*_{0}) in the plasma membrane.

## 3. Results

We first study the dependence of the deformed shape on magnitude of the line force *F*, for the case where force is equal and opposite on the vesicle and plasma membrane. The location of force application is fixed at on both the vesicle and plasma membrane, as shown in figure 3. We vary the strength of the line force in the range of 5–20 in dimensionless terms, which is equivalent to a total force between 66 and 266 pN. In our axisymmetric model, this force represents the total force exerted by all the SNARE complexes attached to the vesicle. For example, Gao *et al.* [11] found that the unzipping force for a single SNARE complex ranges from 2 to 20 pN. This force range should be treated as the lower estimate of the forces exerted *in vivo*, as in this study the SNARE complex unzippering occurs at the range of seconds, whereas the process of neuronal fusion occurs at the scale of microseconds. If we assume that the force of 20 pN or somewhat higher unzips a single SNARE complex *in vivo*, the force range explored in our study would be sufficient to separate three to six SNARE complexes. Figure 3 shows the deformed shapes of the membranes for four different values of . Note that the minimum separation between the membranes does not occur at the bottom of the vesicle, as is usually assumed, but near the point of load application where the SNARE complex is located. For practical purposes, the point of load application can be used as an estimate for the minimum separation. This separation is representative of the separation between the C-termini transmembrane domain of Syb and Syx. Figure 4 shows this estimate of minimum separation versus the applied force. The separation decreases rapidly with increasing applied force, with the rate of decrease of separation becoming slower at higher loads. Note that the minimum separation in the simulations is greater than twice the Debye length of our system, which is consistent with the DH approximation.

As shown in figure 5*a*, the entire vesicle membrane is under tension and the maximum tension occurs at the bottom of the deformed vesicle. However, due to the direction of the applied force, part of the plasma membrane can be under compression despite the pretension. Our results show that the maximum compression occurs at the centre of the plasma membrane as shown in figure 5*b*. Because the applied force in the plasma membrane has a tangential component, the tension is discontinuous across the line of load application as shown in figure 5*b*. Figure 6 shows that these maxima are quite insensitive to the surface charge densities, that is, *increasing the surface charge density on the plasma membrane by a factor approximately 3 does not affect the maximum tension and compression for a given applied force*. The rupture strength of a lipid bilayer is approximately 10 mN m^{−1} [38,39], which on our non-dimensional scale turns out to be approximately 44 units. Therefore, the range of SNARE forces used in our analysis is not sufficient to cause rupture.

The results in figures 3⇑⇑–6 are for membranes that are subjected to equal and opposite forces. We also carried out calculations for the case where the applied forces are always normal to the deformed surfaces. Our numerical result shows that except for the fact that the plasma membrane has much less compression, there are no qualitative differences between these two loading configurations; therefore, plots similar to figures 3⇑⇑–6 for these other loading conditions are given in the electronic supplementary material. For the rest of this paper, we will focus on the case in which the applied forces are equal and opposite.

Next, we vary the location of the SNARE-machinery on both vesicle and plasma membrane () and keep fixed. It is interesting to note that, because the magnitude of the line force is fixed, increasing to is equivalent to increasing the number of SNAREs by a factor of in our simulation. Figure 7 shows that increasing increases the deformation of the membranes, bringing them closer. The tension (compression) in the membranes also increases with , as shown in figures 8 and 9. Note that as increases, there is a smaller jump in tension in the plasma membrane. This result is due to the fact that, as increases, the deformed plasma membrane surface reorients so that the direction of the SNARE force is closer to the surface normal. Figure 10 shows that, as the SNARE-machinery moves away from the centre (increasing ), the maximum tension increases rapidly. Recall in figure 6, the maximum tension and compression vary approximately linearly with the magnitude of the applied force. Therefore, the slope of the lines in figure 10 increases with .

Our pretension boundary condition allows the plasma membrane to deflect. In the real system of synaptic vesicle fusion, pretension in the membranes is due to the osmotic pressure inside the neuron cell relative to its surroundings. In particular, the case of infinite pretension is equivalent to a rigid substrate which was studied earlier by Blount *et al.* [29] (see Introduction).

To determine the influence of pretension in the plasma membrane, we keep the SNARE-machinery force at a fixed magnitude and location, while varying the pretension imposed at the far end of the plasma membrane. Figure 11 shows the deformed shapes for two different values of pretension , the deformation in the plasma membrane as well as that of the vesicle, decreases as the pretension increases. The distributions of tension in the membranes are shown in figure 12 for two different pretensions. For the larger pretension, the entire plasma membrane is under tension. The critical dimensionless pretension where this occurs is 7, as shown in figure 13, where we plot the maximum tension (compression) of the plasma membrane versus pretension. Interestingly, the maximum tension in the vesicle decreases with pretension. Figure 13 shows why this is the case: the minimum separation between the plasma membrane and the vesicle increases as the pretension becomes higher. The increase in separation results in a decrease in the electrostatic repulsion. Because the electrostatic force has a tangential component along the vesicle surface, lowering this force lowers the membrane tension.

To determine the influence of osmotic pressure across the vesicle membrane, we keep the SNARE-machinery force at a fixed magnitude and location, while varying the osmotic pressure across the vesicle membrane. In the synaptic vesicles, the osmolarity is approximately 0.3 Osm, which on non-dimensional scale induces a pressure of about 70 inside the vesicle [40,41]. As stated in [41], the actual osmolarity could be lower than 0.3 Osm, as some of the neurotransmitters can bind with the matrix inside the synaptic vesicle, resulting in a lower osmotic pressure. Figure 14 shows the deformed shapes for three different values of osmotic pressure difference; it is evident that the deformation in the plasma membrane increases as the osmotic pressure increases. However, the minimum separation between the vesicle and the plasma membrane is nearly insensitive to the osmotic pressure. The distributions of tension in the membranes are shown in figure 15 for three different pressure values. As expected, the profile of tension in the vesicle translates upward with a constant value proportional to the osmotic pressure. As a result, the maximum tension increases more rapidly in the vesicle, as compared with the plasma membrane, as shown in figure 16. If the osmotic pressure of the vesicle is kept at a sufficiently high value, so that the peak tension at the vesicle bottom exceeds the failure tension value (approx. 44 units on the dimensionless scale), then it could expose a site for fusion of the membranes.

## 4. Effect of hemifusion

Numerous theoretical and experimental [2,42–45] studies suggested an intermediate step in the fusion process that involves the formation of a stalk and hemifusion diaphragm, followed by membrane rupture. The formation of a stalk implies that instability developed at some location between the two bilayers, which induces a localized lipid rearrangement. Eventually, this rearrangement of lipids leads to the formation of a single bilayer, which is commonly known as hemifusion diaphragm. As mentioned above, our numerical results indicate that the shortest gap between the two membranes lies near the SNARE location. This supports the experimental observation that the transmembrane segment of SNARE induces distortion in the lipid packing around it [46]. It is possible that this distortion will eventually lead to the lipid rearrangement in the membranes resulting in the formation of hemifusion diaphragm. In the following, we propose that, at this point of proximity, hemifusion is initiated and propagates from there towards the bottom of the vesicle.

The exact mechanisms behind this complex sequence of events are still an active area of research. It is difficult to understand hemifusion using continuum membrane theory, because hemifusion is controlled by events on the nanometre scale. In this section, we consider a simpler question: assuming that there is a hemifused state, what is the tension in the hemifused region?

Because hemifusion involves very close contact between the lipid bilayers, it is reasonable to assume that the two membrane surfaces are completely dehydrated. Consistent with this assumption, hemifusion is modelled in our continuum approach by turning off the repulsive force inside a specified region on the membranes, which we define as the hemifused region. Thus, in the hemifused region, the outer surface of the vesicle is free of electrostatic repulsion and hence can be modelled as traction free. A limitation of our model is that it does not account for the fact that the hemifused region has only one bilayer instead of two. Because the real system consists of only one layer and we have neglected concentrated moments and shear at the edge of the hemifused region, the tension predicted by our model is expected to be lower than in the actual situation.

Details of these simulations are presented in the electronic supplementary material. Here, we summarize the main results. Our key result is that the membrane tension increases with the length of the hemifused region. The deformation and tension are more severe than the case of without hemifusion. Our results suggest that hemifusion propagation increases the chances of membrane rupture. Eventually, a fusion pore formation is possible in the hemifused region, where tension is highest.

## 5. Summary and conclusion

In this article, we use a coarse-grained continuum model based on membrane theory to study the deformation and tension in the vesicle and plasma membranes during docking. The zipping action of the SNARE complex is represented by two sets of equal and opposite line forces acting on the membranes. The repulsive interaction between the membranes is represented using DH theory. The magnitude and location of these line forces are varied in our simulations to study their effects on tension and membrane deformation. We also study hemifusion by turning off the repulsive force in a region where the membranes are closest to each other. Our results can be summarized as follows:

— The closest approach between the vesicle and plasma membrane does not occur at the bottom of the vesicle, but near the location where the components of the SNARE complex are inserted into the membranes.

— The maximum in plane tension occurs at the bottom of the vesicle membrane.

— The maximum in plane tension increases linearly with the applied force.

— The maximum in plane tension is insensitive to the charge density of the surfaces. This can be explained by the fact that the change in repulsive forces alters the separation between the membranes, while keeping the local deformation nearly the same as before.

— For small pretensions, the plasma membrane can be under compression, and the maximum compression occurs at the centre, directly below the bottom of the vesicle. For normalized pretension greater than 7, the entire plasma membrane is under tension.

— As the location of force is moved away from the bottom of the vesicle (which is equivalent to increasing the number of SNAREs in our simulation), the tension in the vesicle membrane increases.

— Hemifusion causes an increase in the in-plane tension of both the vesicle and the plasma membrane. We expect that this increase in the magnitude of the tension will eventually lead to rupture of the membranes, leading to fusion pore formation.

## Funding statement

C.Y.H., T.L., P.S., M.B. and A.J. acknowledge the support by the National Institutes of Health under award no. R01 MH099557.

- Received October 11, 2014.
- Accepted November 10, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.