Allosteric stabilization of calcium and phosphoinositide dual binding engages several synaptotagmins in fast exocytosis
Abstract
Synaptic communication relies on the fusion of synaptic vesicles with the plasma membrane, which leads to neurotransmitter release. This exocytosis is triggered by brief and local elevations of intracellular Ca^{2+} with remarkably high sensitivity. How this is molecularly achieved is unknown. While synaptotagmins confer the Ca^{2+} sensitivity of neurotransmitter exocytosis, biochemical measurements reported Ca^{2+} affinities too low to account for synaptic function. However, synaptotagmin’s Ca^{2+} affinity increases upon binding the plasma membrane phospholipid PI(4,5)P_{2} and, vice versa, Ca^{2+} binding increases synaptotagmin’s PI(4,5)P_{2} affinity, indicating a stabilization of the Ca^{2+}/PI(4,5)P_{2} dualbound state. Here, we devise a molecular exocytosis model based on this positive allosteric stabilization and the assumptions that (1.) synaptotagmin Ca^{2+}/PI(4,5)P_{2} dual binding lowers the energy barrier for vesicle fusion and that (2.) the effect of multiple synaptotagmins on the energy barrier is additive. The model, which relies on biochemically measured Ca^{2+}/PI(4,5)P_{2} affinities and protein copy numbers, reproduced the steep Ca^{2+} dependency of neurotransmitter release. Our results indicate that each synaptotagmin engaging in Ca^{2+}/PI(4,5)P_{2} dualbinding lowers the energy barrier for vesicle fusion by ~5 k_{B}T and that allosteric stabilization of this state enables the synchronized engagement of several (typically three) synaptotagmins for fast exocytosis. Furthermore, we show that mutations altering synaptotagmin’s allosteric properties may show dominantnegative effects, even though synaptotagmins act independently on the energy barrier, and that dynamic changes of local PI(4,5)P_{2} (e.g. upon vesicle movement) dramatically impact synaptic responses. We conclude that allosterically stabilized Ca^{2+}/PI(4,5)P_{2} dual binding enables synaptotagmins to exert their coordinated function in neurotransmission.
Editor's evaluation
The calcium dependence of vesicle exocytosis at synapses is a power law with an exponent n = 3 or 4, however, the molecular mechanisms that underpin this highly nonlinear dependence on calcium are unclear. To shed light on this fundamental question the authors build a model where 2 calcium ions bind to the protein synaptotagmin and synaptotagmin binds to the negatively charged lipid PIP2 in the presynaptic membrane. Simulations fit best the data from the calyx of Held synapse when 3 synaptotagmin molecules each bind calcium and PIP2. This compelling model shows that each CasynaptotagminPIP2 complex reduces the energy barrier for vesicle fusion by ~5k, thus, fast exocytosis at CNS synapses may require only 3 CasynaptogaminPIP2 molecules to achieve submillisecond speeds of vesicle fusion.
https://doi.org/10.7554/eLife.74810.sa0eLife digest
For our brains and nervous systems to work properly, the nerve cells within them must be able to ‘talk’ to each other. They do this by releasing chemical signals called neurotransmitters which other cells can detect and respond to.
Neurotransmitters are packaged in tiny membranebound spheres called vesicles. When a cell of the nervous system needs to send a signal to its neighbours, the vesicles fuse with the outer membrane of the cell, discharging their chemical contents for other cells to detect. The initial trigger for neurotransmitter release is a short, fast increase in the amount of calcium ions inside the signalling cell. One of the main proteins that helps regulate this process is synaptotagmin which binds to calcium and gives vesicles the signal to start unloading their chemicals.
Despite acting as a calcium sensor, synaptotagmin actually has a very low affinity for calcium ions by itself, meaning that it would not be efficient for the protein to respond alone. Synpatotagmin is more likely to bind to calcium if it is attached to a molecule called PIP_{2}, which is found in the membranes of cells The effect also occurs in reverse, as the binding of calcium to synaptotagmin increases the protein’s affinity for PIP_{2}. However, how these three molecules – synaptotagmin, PIP_{2}, and calcium – work together to achieve the physiological release of neurotransmitters is poorly understood.
To help answer this question, Kobbersmed, Berns et al. set up a computer simulation of ‘virtual vesicles’ using available experimental data on synaptotagmin’s affinity with calcium and PIP_{2}. In this simulation, synaptotagmin could only trigger the release of neurotransmitters when bound to both calcium and PIP_{2}. The model also showed that each ‘complex’ of synaptotagmin/calcium/PIP_{2} made the vesicles more likely to fuse with the outer membrane of the cell – to the extent that only a handful of synaptotagmin molecules were needed to start neurotransmitter release from a single vesicle.
These results shed new light on a biological process central to the way nerve cells communicate with each other. In the future, Kobbersmed, Berns et al. hope that this insight will help us to understand the cause of diseases where communication in the nervous system is impaired.
Introduction
Regulated neurotransmitter (NT) release from presynaptic terminals is crucial for information transfer across chemical synapses. NT release is triggered by action potentials (APs), which are transient de and repolarizations of the presynaptic membrane potential that induce Ca^{2+} influx through voltagegated channels. The resulting brief and local elevations of the intracellular Ca^{2+} concentration ([Ca^{2+}]_{i}) trigger the fusion of NTcontaining synaptic vesicles (SVs) from the socalled readily releasable pool (RRP), whose SVs are localized (docked) at the plasma membrane and molecularly matured (primed) for fusion (Kaeser and Regehr, 2017; Südhof, 2013; Verhage and Sørensen, 2008). A high Ca^{2+} sensitivity of NT release is needed to achieve fast responses to the very short APinduced Ca^{2+} transient and correspondingly, the SV fusion rate depends to the 4th5th power on the [Ca^{2+}]_{i} (Bollmann et al., 2000; Burgalossi et al., 2010; Heidelberger et al., 1994; Schneggenburger and Neher, 2000). Accordingly, previous models of NT release have assumed the successive binding of five Ca^{2+} ions to a sensor that regulates release (Bollmann et al., 2000; Lou et al., 2005; Schneggenburger and Neher, 2000). However, how these macroscopic properties arise from the molecular components involved in SV fusion is still unknown.
The energy for SV fusion is provided by the assembly of the neuronal SNARE complex, which consists of vesicular synaptobrevin/VAMP and plasma membrane bound SNAP25 and syntaxin proteins (Jahn and Fasshauer, 2012; Südhof, 2013). Ca^{2+} sensitivity of SV fusion is conferred by the vesicular protein synaptotagmin (syt), which interacts with the SNAREs (Brewer et al., 2015; Littleton et al., 1993; Mohrmann et al., 2013; Schupp et al., 2017; Zhou et al., 2015; Zhou et al., 2017). Several syt isoforms are expressed in synapses. Depending on the synapse type (e.g. mouse hippocampal pyramidal neurons or the Calyx of Held), syt1 or syt2 is required for synchronous, Ca^{2+}induced fusion (Geppert et al., 1994; Kochubey et al., 2016; Kochubey and Schneggenburger, 2011; Südhof, 2013). These two syt isoforms are highly homologous and contain two cytosolic, Ca^{2+}binding domains, C2A and C2B (Südhof, 2002), of which the C2B domain has been shown to be essential, and in some cases even sufficient, for synchronous NT release (Bacaj et al., 2013; Gruget et al., 2020; Kochubey and Schneggenburger, 2011; Lee et al., 2013; Mackler et al., 2002). The C2B domain contains two Ca^{2+} binding sites on its top loops (Fernandez et al., 2001). In addition, a second binding site allows the C2B domain to bind to the signaling lipid phosphatidylinositol 4,5bisphosphate (PI(4,5)P_{2}) located in the plasma membrane (Bai et al., 2004; FernándezChacón et al., 2001; Honigmann et al., 2013; Li et al., 2006; Xue et al., 2008), but might also participate in (possibly transient) SNARE interactions (Brewer et al., 2015; Zhou et al., 2015; Zhou et al., 2017). A third site, located in the far end of the C2B domain (R398 and R399 in mouse syt1), is also involved in both SNARE and membrane contacts (Nyenhuis et al., 2021; Xue et al., 2008; Zhou et al., 2015). Via these interactions, the syt C2B domain can induce close membranemembrane contact in vitro (Araç et al., 2006; Chang et al., 2018; Honigmann et al., 2013; Nyenhuis et al., 2021; Seven et al., 2013; Xue et al., 2008), stable vesiclemembrane docking (Chang et al., 2018; Chen et al., 2021; de Wit et al., 2009), as well as dynamic vesiclemembrane association upon Ca^{2+} influx into the cell (Chang et al., 2018).
Despite its central role as the Ca^{2+} sensor for NT release, the intrinsic Ca^{2+} affinity of the isolated syt C2B domain is remarkably low (KD ≈ 200 µM, Radhakrishnan et al., 2009; van den Bogaart et al., 2012), much lower than the Ca^{2+} sensitivity of NT release (Bollmann et al., 2000; Schneggenburger and Neher, 2000). However, binding of the C2B domain to PI(4,5)P_{2}, which is enriched at synapses (van den Bogaart et al., 2011a), drastically increases its Ca^{2+} affinity (van den Bogaart et al., 2012). Similarly, the affinity for PI(4,5)P_{2} increases upon Ca^{2+} binding (PérezLara et al., 2016; van den Bogaart et al., 2012). This indicates a positive allosteric coupling between the binding sites for Ca^{2+} and PI(4,5)P_{2}, which promotes dual binding of Ca^{2+}/PI(4,5)P_{2} (Li et al., 2006; Radhakrishnan et al., 2009; van den Bogaart et al., 2012). As binding of both molecules to syt is involved in fusion (Kedar et al., 2015; Li et al., 2006; Mackler et al., 2002; Mackler and Reist, 2001; Wang et al., 2016; Wu et al., 2021a; Wu et al., 2021b), this positive allosteric coupling might be central to syt’s function in triggering Ca^{2+}induced exocytosis (van den Bogaart et al., 2012).
In this paper, we developed a mathematical model in which the dual binding of the C2B domain to Ca^{2+} and PI(4,5)P_{2} promotes fusion. The model, which is based on the measured affinities and allostericity of Ca^{2+} and PI(4,5)P_{2} binding, describes stochastic binding/unbinding reactions at the level of individual syts and stochastic SV fusion events. The model predicts that each C2B domain engaging in dual Ca^{2+}/PI(4,5)P_{2} binding lowers the energy barrier for fusion by ~5 k_{B}T. Our results indicate that during fast NT release most fusion events occur once three syts per SV simultaneously engage their C2B domains in dual Ca^{2+}/PI(4,5)P_{2} binding. This simultaneous engagement of multiple syts crucially relies on the positive allosteric coupling between Ca^{2+} and PI(4,5)P_{2} binding. We explored consequences of putative mutations affecting Ca^{2+}/PI(4,5)P_{2} binding and/or the allosteric coupling between the binding sites of both species and suggest that changes of the allostericity contribute to dominant negative effects. Moreover, dynamic changes of PI(4,5)P_{2} accessibility for the syts (e.g. induced by SV movement to the plasma membrane) are predicted to dramatically impact synaptic responses. We conclude that allosterically stabilized Ca^{2+}/PI(4,5)P_{2} dual binding to the C2B domain forms the molecular basis for synaptotagmins to exert their cooperative control of neurotransmitter release.
Results
An experimentbased model of the triggering mechanism for SV fusion based on molecular interactions between syt, Ca^{2+}, and PI(4,5)P_{2}
To develop an experimentbased model of NT release based on molecular properties of syt, we first described the reaction scheme of a single C2B domain. The C2B domain binds PI(4,5)P_{2} and two Ca^{2+} ions (FernándezChacón et al., 2002; Honigmann et al., 2013; Mackler et al., 2002; van den Bogaart et al., 2012; Xue et al., 2008). We assumed the simplest case of simultaneous association of both Ca^{2+} ions to the syt1 C2B domain. Therefore, in our model the C2B domain can be in four different states (Figure 1A): (1) an ‘unbound’ state, (2) a PI(4,5)P_{2}bound state, (3) a state with two Ca^{2+} ions bound, and (4) a ‘dualbound’ state in which the C2B simultaneously engages Ca^{2+}/PI(4,5)P_{2} binding. The affinities for Ca^{2+} and PI(4,5)P_{2} were set to those measured in vitro (K_{D,2Ca2+} and K_{D,PIP2})(van den Bogaart et al., 2012). Binding of PI(4,5)P_{2} to the C2B domain was shown to increase the domain’s affinity for Ca^{2+} and vice versa, indicating a positive allosteric coupling between the two binding sites (van den Bogaart et al., 2012). We therefore implemented a positive allosteric stabilization of the dualbound state in the model (illustrated by the red shaded areas of the C2B domain in Figure 1A) by introducing the allosteric factor (A=0.00022, see Materials and methods van den Bogaart et al., 2012) which slows down the Ca^{2+} and PI(4,5)P_{2} dissociation from the dualbound state.
We next extended the model to the level of the complete SV. On average, SVs contain 15 copies of syt1 (Figure 1B; Takamori et al., 2006), which may work together to regulate SV fusion. Spontaneous release occurs at low rates (with a rate constant ‘L_{+}’), reflected by a high initial energy barrier for SV fusion (Figure 1C). Because syt’s stimulation of SV fusion likely relies on the simultaneous binding of both Ca^{2+} and PI(4,5)P_{2} (Kedar et al., 2015; Li et al., 2006; Mackler et al., 2002; Mackler and Reist, 2001; Wang et al., 2016; Wu et al., 2021a; Wu et al., 2021b), we assumed that each dualbound C2B domain promotes exocytosis by lowering this barrier. How this might be achieved exactly is unknown, but could involve bridging plasma and SV membranes (Figure 1A), changing the curvature of the plasma membrane (Figure 1—figure supplement 1), changing the local electrostatic environment, or directly or indirectly promoting SNARE complex assembly (Bhalla et al., 2006; Martens et al., 2007; Ruiter et al., 2019; Schupp et al., 2017; Tang et al., 2006; van den Bogaart et al., 2011b; Zhou et al., 2015; Zhou et al., 2017). We assumed that multiple syts may progressively lower the energy barrier by the successive engagement of their C2B domains in dual Ca^{2+}/PI(4,5)P_{2} binding, and investigated the simplest scenario, in which each dualbound C2B domain contributed the same amount of energy (E_{syt}) thereby increasing the fusion rate by the same factor (f) (Figure 1C; Schotten et al., 2015). The number of syts that can simultaneously promote fusion may be limited by their access to PI(4,5)P_{2} in the plasma membrane. This could be due to limited space beneath the SV or limited molecular access to bind PI(4,5)P_{2} clusters and/or SNAREs (de Wit et al., 2009; Mohrmann et al., 2013; Rickman and Davletov, 2003). We therefore created a model that stochastically describes the binding status of individual SVs (Figure 1—figure supplement 2) with a limited number of PI(4,5)P_{2} association possibilities (‘slots’) and investigated how this number affects physiological responses.
At least three PI(4,5)P_{2} binding slots are required to reproduce release kinetics from the calyx of Held synapse
A hallmark of the NT release reaction is its large dynamic range in response to Ca^{2+} stimuli as impressively demonstrated by experimental data from the calyx of Held synapse where release latencies (defined as the time of the fifth SV fusion after the stimulus) and exocytosis rates have been measured for a broad range of Ca^{2+} concentrations using Ca^{2+} uncaging (Bollmann et al., 2000; Kochubey and Schneggenburger, 2011; Lou et al., 2005; Schneggenburger and Neher, 2000; Sun et al., 2007). At this wellestablished model synapse, fast NT release is controlled by syt2, which is functionally redundant with syt1 in neurons (Kochubey et al., 2016; Xu et al., 2007).
We evaluated whether our model could reproduce this Ca^{2+} dependence by simulating release latencies and peak release rates in response to steplike Ca^{2+} stimuli. The ability to reproduce the experimental data depended on the number of ‘slots’ for syt PI(4,5)P_{2} binding (Figure 2A). We first fitted the free parameters in our model by optimizing the agreement (i.e. by reducing a predefined cost function, see Materials and methods) between model predictions and release rates and latencies determined experimentally by Kochubey and Schneggenburger (Figure 2A; Kochubey and Schneggenburger, 2011). During this fitting process, we took the entire distribution of the experimentally obtained release latencies into account by using the likelihood function (see Materials and methods). This was not feasible for the experimental peak release rates (since accurate computation of the maximum rate of stochastic events is not feasible), which were therefore compared to the closed form solution of the model (see Materials and methods). Because the affinities for Ca^{2+} and PI(4,5)P_{2}, and the allosteric coupling between both species (K_{D,2Ca2+}, K_{D,PIP2} and A, Figure 1A, Figure 1—figure supplement 2) as well as the RRP size (including its variance) were taken from literature (Table 1; Figure 2—figure supplement 1; van den Bogaart et al., 2012; Wölfel and Schneggenburger, 2003), we only had to estimate five parameters: (1) the binding rate constant of the two Ca^{2+} ions (α), (2) the binding rate constant of PI(4,5)P_{2} (γ), (3) the PI(4,5)P_{2} concentration ([PI(4,5)P_{2}]), (4) the reduction of the energy barrier for fusion induced by Ca^{2+}/PI(4,5)P_{2} dual binding of one C2B domain (${E}_{syt}$) and (5) a fixed delay (d) between time of uncaging and response onset, like in previous models (see Kochubey and Schneggenburger, 2011; Schneggenburger and Neher, 2000). Having obtained the best fit parameters, peak release rates and release latencies (Figure 2A) were stochastically simulated based on the closedform solution of the Markov model (see Materials and methods).
We systematically varied the number of slots (M_{slots}) from one to six, optimized the free parameters for each of these choices, and compared the best fit solutions. With a single slot (M_{slots} = 1), the model accounted for experimentally observed release latencies, but failed to reproduce the steep relationship between [Ca^{2+}]_{i} and peak release rates (Figure 2A). Adding more slots strongly improved the agreement with experimental data. The best agreement was found with three slots (M_{slots} = 3) and the agreement slightly decreased with more slots (M_{slots} >3) (Figure 2B). However, all models with at least three slots (M_{slots} ≥3) captured the steep dependency of peak release rates on [Ca^{2+}]_{i}, with a maximum slope of 4–5 on a doublelogarithmic plot (Figure 2C; Schneggenburger and Neher, 2000).
Our model made it possible to inspect the fate of each individual fusing SV, including the number of synaptotagmins dually binding Ca^{2+} and PI(4,5)P_{2} just before fusion. Remarkably, even in models with more than three slots (M_{slots} >3), fast NT release, induced by moderate to high [Ca^{2+}]_{i}, was predicted to primarily engage three dualbound syt C2B domains. At lower [Ca^{2+}]_{i} fewer dualbound syts mediated fusion (Figure 2D–E, Figure 2—figure supplement 2, Figure 2—figure supplement 3). Correspondingly, the estimated reduction of the energy barrier for fusion by each Ca^{2+}/PI(4,5)P_{2} dual binding C2B domain was similar (E_{syt}≈5 k_{B}T) for all versions of the model with at least three slots (M_{slots} ≥3) (Figure 2F, see Table 2 for best fit model parameters for each setting of M_{slots}). Our model thus predicts that most fusion events occur when three syts actively reduce the energy barrier for fast SV fusion.
We next investigated to which extent the estimated number of syts working together for fusion depended on some of the assumptions of our model. For instance, if each syt dualbinding Ca^{2+}/PI(4,5)P_{2} had a lower effect on the vesicle fusion rate, might this be compensated by more syts working together during fusion? We investigated this by manually forcing a lower individual contribution to the energy barrier for fusion in a model with six slots and refitting the remaining parameters. Under such conditions, the dependence of the peak release rate on [Ca^{2+}]_{i} became too steep, indicating that too many syts working together make the Ca^{2+} sensitivity unnaturally high (Figure 2—figure supplement 4). We then investigated how the assumption of simultaneous binding of two Ca^{2+} ions to the C2B domain affected the conclusions. If each C2B domain only bound a single Ca^{2+} ion, the dependence of the peak release rate on [Ca^{2+}]_{i} was too shallow, even in a model with six slots. Allowing the number of Ca^{2+} ions binding to one C2B domain to vary in a macroscopic version of the model with six slots predicted the binding of 1.53 Ca^{2+} ions per C2B domain and most NT release commencing with four or fewer cooperating syts (Figure 2—figure supplement 5). This confirms that most C2B domains need to bind two Ca^{2+} ions to exert their effect. Simulating a model with consecutive Ca^{2+} binding to the two binding sites of the C2B domain on all syts of RRP SVs would be computationally too costly (and would involve additional, unknown parameters). However, we show that our simplification of simultaneous binding aligns with such a more complex model if the binding of the second Ca^{2+} ion is favored (Figure 2—figure supplement 6), which is reasonable based on the proximity to negatively charged lipid headgroups following the insertion of the C2B domain into the plasma membrane. Thus, while our model is a simplification of the reality, the main conclusion on the number of slots needed (M_{slots} ≥3) and the number of syts sufficient to mediate fast NT release (≤4), are robust estimates. Because our molecular model assuming three slots (M_{slots} = 3) and simultaneous binding of two Ca^{2+} ions to the syt C2B domains was the simplest model that accounted for the experimental data, we used it for further simulations.
The best fit parameters for three slots revealed rapid association rate constants for Ca^{2+} and PI(4,5)P_{2} to the C2B domain and PI(4,5)P_{2} levels corresponding to a concentration of ~1 µM in an in vitro setting (Table 2). Predicted responses obtained using the best fit parameters were sensitive to changes of either of these parameters (Figure 2—figure supplement 7). For instance, higher levels of PI(4,5)P_{2} decreased the release latencies and increased the rate of fusion, and changing the Ca^{2+} association rate constant (α) affected the release latencies much more than changing the PI(4,5)P_{2} association rate constant (γ). We verified that these parameters represent unique solutions by systematically exploring the parameter space with f (which relates to the lowering of the fusion barrier for each syt dualbinding Ca^{2+}/PI(4,5)P_{2}) fixed to the best fit value (f=128), which revealed a clear minimum at the best fit parameters (Figure 2G, darkest ball). We furthermore confirmed that this f value was optimal by systematically varying f and fitting all other parameters (Figure 2H).
The number of syt proteins preassociated to PI(4,5)P_{2} at rest influences the SV’s Ca^{2+} responsiveness
The steady state concentration of PI(4,5)P_{2} determines the probability of syts associating to PI(4,5)P_{2} at rest. With the best fit parameters, our model predicts that at rest ([Ca^{2+}]_{i}=50 nM) most SVs associate to PI(4,5)P_{2} by engaging one (~42%), two (~33%) or three (~8%) syts (Figure 3A, see Figure 3—figure supplement 1 for behavior in the model with M_{slots} = 6). With a steplike Ca^{2+} stimulus to 50 µM, SVs with two or three preassociated syts mediated most of the fastest (<0.5ms) SV fusions (Figure 3B). Consequently, changing the steady state PI(4,5)P_{2} concentration (which changes the number of preassociated syts/SVs) largely impacted the release latencies (defined as the timing of the fifth SV that fuses, Kochubey and Schneggenburger, 2011; Figure 2—figure supplement 7). Due to the allosteric coupling, the Ca^{2+} affinities of syts prebound to PI(4,5)P_{2} are increased and SVs with more PI(4,5)P_{2} interactions are more responsive to the Ca^{2+} stimulus (Figure 3C). Thus, at the single SV level, the number of preassociated syts to PI(4,5)P_{2} at rest plays a role in very fast (submillisecond) SV release and causes heterogeneity in release probability among RRP SVs.
Allosteric stabilization of Ca^{2+}/PI(4,5)P_{2} dual binding is necessary to synchronize multiple C2B domains for fast SV fusion
An important feature of our model is the inclusion of a positive allosteric interaction between Ca^{2+} and PI(4,5)P_{2} binding to the C2B domain which we based on increased affinities measured in vitro (van den Bogaart et al., 2012). To explore the physiological relevance of this allostericity, we investigated how individual SVs engaged their syt C2B domains in Ca^{2+}/PI(4,5)P_{2} dual binding in response to a stepwise increase of [Ca^{2+}]_{i} to 50 µM (Figure 4A) with or without this allosteric coupling. We did this by following the fate of the RRP SVs in stochastic simulations, four of which are illustrated in Figure 4B. Under normal conditions (with allostericity), syt C2B domains quickly associated both Ca^{2+} and PI(4,5)P_{2} and their respective allosteric stabilization slowed the dissociation of either species resulting in a lifetime of their dual binding of ~1.3 ms on average. This enabled the successive engagement of three dualbound C2B domains for most RRP SVs (including all four illustrated SVs). The average waiting time for fusion for the RRP SVs was ~1.1 ms (Figure 4B, fusion indicated by circles). By inspecting the average behavior of the entire RRP of SVs it became clear that the overall release rate closely followed the population of SVs engaging three C2B domains in dual Ca^{2+}/PI(4,5)P_{2} binding, illustrating the importance of engaging three syts for fast SV fusion in this model (Figure 4C). We also simulated the postsynaptic response produced by this NT release by convolving the SV release rate with a typical postsynaptic response to the fusion of a single SV (see Materials and methods), which revealed synchronous and large Excitatory Post Synaptic Currents (EPSCs)(Figure 4D).
We then explored what would happen without the allosteric stabilization of Ca^{2+}/PI(4,5)P_{2} dual binding (by setting A=1; Figure 4E). In this case, the C2B domains still quickly associated Ca^{2+} and PI(4,5)P_{2}, but without the allosteric slowing of Ca^{2+}/PI(4,5)P_{2} dissociation the lifetime of dualbound C2B domains was dramatically reduced to an average of ~0.0003 ms. This made it very improbable to engage multiple C2B domains in dual Ca^{2+}/PI(4,5)P_{2} binding (Figure 4E). In turn, without the simultaneous engagement of multiple syts dualbinding Ca^{2+}/PI(4,5)P_{2}, NT release became very unlikely. In fact, none of the randomly chosen four RRP SVs fused within 4 ms (Figure 4E). Inspection of the average behavior of the entire RRP revealed that only few SVs engaged more than one syt C2B domain in dual Ca^{2+}/PI(4,5)P_{2} binding, resulting in a very low fusion rate (Figure 4F). Correspondingly, postsynaptic EPSCs were severely disrupted, and most release events were illsynchronized single SV fusion events (Figure 4G). It was furthermore not possible to fit a model without the allosteric stabilization to the experimental dataset (Figure 4—figure supplement 1). Thus, the positive allosteric coupling between Ca^{2+} and PI(4,5)P_{2} is fundamental for the syts to simultaneously and persistently engage multiple C2B domains per SV in Ca^{2+}/PI(4,5)P_{2} dual binding.
Many syts per SV speed up exocytosis by increasing the probability of Ca^{2+}/PI(4,5)P_{2} dual binding
Our model suggests that only a few syts simultaneously binding Ca^{2+} and PI(4,5)P_{2} are required to promote fast SV fusion (Figure 2). Yet, a total of 15 copies are expressed per SV on average (Takamori et al., 2006), which raises the question why SVs carry such excess and whether and how the additional syt copies contribute to the characteristics of Ca^{2+}induced synaptic transmission. To investigate this, we simulated Ca^{2+} uncaging experiments with reduced numbers of syts per SV while keeping all other parameters in the model constant. With fewer syts, release latencies increased and peak release rates reduced. Defects were particularly prominent for reductions to less than three copies per SV (Figure 5A). Further exploration indicated that the responses slowed down upon reductions in syt copy number because it took SVs longer to simultaneously engage three C2B domains in dual Ca^{2+}/PI(4,5)P_{2} binding and that fewer SVs reached this state (Figure 5B).
While Ca^{2+} uncaging stimuli are exquisitely suited to map the full range of synaptic responses, synaptic transmission is physiologically triggered by APs that induce shortlived Ca^{2+} transients. To stochastically predict responses to such timevarying Ca^{2+} stimuli, we implemented our model using the Gillespie algorithm. After verifying that this model implementation agreed with the initial implementation (Figure 5—figure supplement 1), we simulated responses to a typical APinduced Ca^{2+} wave that RRP SVs experience (Figure 5C, top panel)(Wang et al., 2008). With 15 syts per SV, APinduced EPSCs were large and synchronous, but reducing their number decreased response amplitudes (Figure 5C). Removal of one syt already reduced the average EPSC amplitude by ~10% and removal of half (7/15) of its copies reduced it by ~72% (Figure 5—figure supplement 2, for representative example traces see Figure 5C). Note, however, that our model only describes the functioning of syt1 /syt2 and therefore does not include other Ca^{2+} sensors, like syt7 and Doc2B, which may mediate release in case of syt1 /syt2 loss (Bacaj et al., 2013; Kochubey et al., 2016; Kochubey and Schneggenburger, 2011; Luo and Südhof, 2017; Sun et al., 2007; Wen et al., 2010; Yao et al., 2011).
As the number of syts per SV has a large impact on fusion kinetics, we wondered to what extent fluctuations in the number of syts per SV affected the variance in APevoked responses in case of their imperfect sorting. Strikingly, however, varying the number of syts per SV over a large range (Poisson distribution with mean = 15, Figure 5D1, E) did not increase the variability of APevoked in synaptic responses while fluctuations of the RRP size strongly impacted them (Figure 5D2, E). This shows that although release kinetics strongly depend on the average number of syts per SV, the system is rather insensitive to fluctuations around this number between individual SVs. Taken together, our data show that although only a subset of syts are required to simultaneously bind Ca^{2+} and PI(4,5)P_{2} to induce fusion, all SV syts contribute to the high rates of NT release by increasing the probability that several syts simultaneously engage in dual Ca^{2+}/PI(4,5)P_{2} binding.
Besides the number of syts, the PI(4,5)P_{2} levels also determine how likely it is for syts to engage in dual Ca^{2+}/PI(4,5)P_{2} binding at an SV (see Figure 3 and Figure 2—figure supplement 7). We therefore reasoned that upregulation of PI(4,5)P_{2} levels, which are dynamically regulated (Jensen et al., 2022), could potentially compensate for reduced syt expression. To investigate this, we refitted the models with reduced syt levels to the experimental Ca^{2+} uncaging data (Kochubey and Schneggenburger, 2011) and only allowed the PI(4,5)P_{2} concentration in the slots ([PI(4,5)P_{2}]) to vary. Strikingly, increasing [PI(4,5)P_{2}] fully rescued the characteristics of NT release upon reductions in syt levels down to 3 syts per SV (corresponding to an 80% reduction) by restoring the number and speed of C2B domains engaging in Ca^{2+}/PI(4,5)P_{2} dual binding (Figure 5—figure supplement 2A–C). The required increase in [PI(4,5)P_{2}] ranged froma factor ~1.1 (14 syts) toa factor ~10 (3 syts, Figure 5—figure supplement 2D). These elevations also fully restored simulated APevoked responses when at least three syts per SV were present (Figure 5—figure supplement 2E, F). Altogether, these data indicate that upregulating [PI(4,5)P_{2}] is a potential, powerful compensatory mechanism to rescue reductions of NT release in case the number of (functional) syts per SV is reduced to no less than three. We note that this compensatory mechanism may strongly influence experimentally observed effects of stoichiometric changes.
Evaluation of mutants affecting Ca^{2+} binding to the C2B domain reveals diverse effects on APevoked transmission
Ca^{2+} sensing of syts depends on negatively charged aspartate (D) sidechains of the C2B domain whose positions are optimal to bind Ca^{2+} ions (Fernandez et al., 2001). The local negative charges of the Ca^{2+} binding sites are reduced/neutralized upon Ca^{2+} binding. The Ca^{2+} binding pockets of the C2B domain have been extensively studied using various mutations (Bradberry et al., 2020; Guan et al., 2017; Kochubey and Schneggenburger, 2011; Lee et al., 2013; Mackler et al., 2002; Nishiki and Augustine, 2004; Shin et al., 2009). Mutations that remove or invert the negative charge of the Ca^{2+} binding sites (by mutation to asparagine (N) or lysine (K), ‘DN’ or ‘DK’) block Ca^{2+} binding and severely reduce exocytosis, even when coexpressed together with the wildtype protein (Bradberry et al., 2020; Kochubey and Schneggenburger, 2011; Lee et al., 2013; Mackler et al., 2002). Other mutations also interfere with Ca^{2+} binding and exocytosis but hold the same pocket charge (e.g. mutation to Glutamate, ‘DE’) (Bradberry et al., 2020). While both types of mutations may similarly interfere with Ca^{2+} binding, they may differentially affect the allosteric mechanism. The allosteric coupling between the Ca^{2+} and PI(4,5)P_{2} binding sites might be (in part) mediated by electrostatic interactions (van den Bogaart et al., 2012), which would imply that the negatively charged Ca^{2+} binding pocket repels PI(4,5)P_{2} until Ca^{2+} reverses the electrostatic charge, and vice versa. Following this assumption, chargealtering mutations within the Ca^{2+} binding pockets (‘DN’, ‘DK’) would partially activate the allosteric coupling mechanism and thereby affect the domain’s PI(4,5)P_{2} affinity (which would not be the case in mutants conserving the charges (‘DE’)). We explored this possibility in our model using two different hypothetical Ca^{2+} site mutants (Figure 6A). We investigated the effect of these mutants on APinduced synaptic transmission under homozygous and heterozygous expression conditions (combined expression of mutant and WT with a total of 15 syts per SV; Figure 6B).
The first mutant, the ‘Ca^{2+}binding’ mutant, had a lower Ca^{2+} affinity (10xK_{D,2Ca2+}), but all other properties were the same as in the WT C2B domain. This might be similar to a mutant with a mutation of the binding pocket which conserves its charge (e.g. ‘DE’). When homozygously expressed, this mutant showed eEPSCs with a~50% reduced amplitude and faster kinetics compared to the WT condition (Figure 6C_{1}, Figure 6D top). Heterozygous expression, with 8 WT and 7 mutant syts per SV, only caused a small decrease in mean eEPSC amplitude compared to the expression of 15 WT syts per SV (Figure 6C1–D, bottom), showing that this mutant is relatively mild.
The second hypothetical mutation was designed to not only abolish Ca^{2+} binding, but to also mimic the Ca^{2+}bound state. Thereby, this mutant featured a high PI(4,5)P_{2} affinity as if the allosteric interaction between Ca^{2+} and PI(4,5)P_{2} was permanently ‘on’. This might represent an extreme example of a mutation electrostatically reducing/inverting the negative charges of the Ca^{2+} binding pocket (e.g. ‘DN’, ‘DK’). We termed this mutant the ‘no Ca^{2+} binding, Aon’ mutant (Figure 6A). This mutant showed no NT release in response to the Ca^{2+} transient in a homozygous condition (Figure 6C2–D, top), which is explained by its inability to bind Ca^{2+}. A major detrimental effect of the mutant was observed when coexpressed with the wildtype protein: When half of the syts on the SV were mutated (heterozygote), the amplitude of simulated eEPSCs was strongly reduced (Figure 6C2–D, bottom). Merely four mutant proteins expressed together with 11 WT proteins already decreased eEPSC amplitudes by ~70% (Figure 6D, bottom), indicating a strong dominant negative effect. The strong inhibition is a result of the mutant’s increased PI(4,5)P_{2} affinity leading to occupation of PI(4,5)P_{2} binding slots on the membrane with this Ca^{2+}insensitive mutant which blocks the association of the Ca^{2+} sensitive and SV fusionpromoting WT proteins. In comparison, a mutant not able to bind Ca^{2+} but having a normal PI(4,5)P_{2} affinity (‘no Ca^{2+} binding, Aoff’ mutant, which could represent a more extreme form of the ‘DE’ mutant) had a much weaker effect (Figure 6—figure supplement 1). This indicates that the allosteric interaction between Ca^{2+} and PI(4,5)P_{2} can play a prominent role in the severity of mutations.
Rapid changes of accessible PI(4,5)P_{2} dramatically impact synaptic shortterm plasticity
In our model, we describe the PI(4,5)P_{2} levels in concentration units, because our model is based on syt affinities for Ca^{2+} and PI(4,5)P_{2} determined in vitro (van den Bogaart et al., 2012). The estimated concentration of PI(4,5)P_{2} not only depends on the local density of PI(4,5)P_{2} in the membrane, but also on the accessibility syt has to PI(4,5)P_{2}. While all species (Ca^{2+}, PI(4,5)P_{2}, and syt C2B) are homogenously accessible in the aqueous solution of the in vitro setting (van den Bogaart et al., 2012), at the synapse the syt C2 domains have constrained motility due to their vesicular association and PI(4,5)P_{2} is restricted to (clusters on) the plasma membrane (Milosevic et al., 2005; van den Bogaart et al., 2011a). This implies that the positioning of SVs with respect to the plasma membrane has an impact on the PI(4,5)P_{2} concentration accessible to syts. We so far assumed that all syts of RRP SVs are exposed to the same PI(4,5)P_{2} levels. This could be the case if all SVs are similarly docked to the plasma membrane. However, when considering more complex stimulation paradigms such as repetitive AP stimulation, this may no longer be valid as several studies reported activitydependent changes in SV positioning on a millisecond timescale (Chang et al., 2018; Kusick et al., 2020; Miki et al., 2016). Rapid changes in accessible PI(4,5)P_{2} may thus contribute to shortterm plasticity, the alteration of responses on a millisecond timescale (Abbott and Regehr, 2004; Kobbersmed et al., 2020; Neher and Brose, 2018; Silva et al., 2021). Recent studies reported that mutations of positively charged amino acids of the C2B domain (lysines, ’Ks’, implicated in binding of PI(4,5)P_{2} and/or the SNAREs and arginines, ‘Rs’, implicated in binding the plasma membrane and/or the SNAREs) resulted in a loss of SV docking and severely reduced neurotransmission (Chang et al., 2018; Chen et al., 2021; Li et al., 2006; Xue et al., 2008). Strikingly, SV docking in these mutants was rapidly restored within milliseconds after an AP which also led to enhanced synaptic transmission in response to a second AP given 10ms after the first (Chang et al., 2018). We explored such a situation in the context of our model by driving exocytosis with two successive APinduced Ca^{2+} transients and assuming either constant PI(4,5)P_{2} levels for syts in wildtype synapses (i.e. all RRP SVs similarly docked) or initially reduced and activitydependent increasing PI(4,5)P_{2} levels for syts in mutant synapses (where SVs docked after the first AP; Figure 7A). We studied the consequence of a mutation that would only affect SV docking at steady state (as may be the case upon mutation of the arginines 398 and 399 of mouse syt 1, ‘R398,399Q’) in (Figure 7). This resulted in a markedly decreased initial response (Figure 7B and C), but repeated activation induced a large facilitation of responses (indicated by a large paired pulse ratio: quotient of the second EPSC amplitudes divided by the first) (Figure 7D). We conclude that dynamic changes in the PI(4,5)P_{2} levels accessible to syts – which may be caused by activitydependent SV relocation – strongly impact synaptic shortterm plasticity. Mutations of the C2B domain that reduce its PI(4,5)P_{2} affinity (as is likely the case upon mutation of the lysine residues 325 and 327 in syt1 or 327, 328 and 332 in syt2) may be more detrimental because even when the effective PI(4,5)P_{2} concentration accessible to syts is restored upon activitydependent SV redocking, syt association to PI(4,5)P_{2} will still be less probable.
Discussion
Here, we propose a quantitative, experimentbased model describing the function of syt in SV fusion on a molecular level based on biochemical properties determined in vitro. In our model, syt acts by lowering the energy barrier for SV fusion by dual binding to Ca^{2+} and PI(4,5)P_{2}. When allowing at least three dualbound syts per SV at a time, this model can explain the steep Ca^{2+} dependence of NT release observed at the calyx of Held synapse (Kochubey and Schneggenburger, 2011). Exploring this model led to the following conclusions:
The positive allosteric interaction between Ca^{2+} and PI(4,5)P_{2} is crucial for fast SV fusion as it stabilizes the dualbound state which allows multiple syts to successively lower the energy barrier for SV fusion;
At least three slots per SV for syt Ca^{2+}/PI(4,5)P_{2} dual binding are needed to achieve the speed and Ca^{2+} sensitivity inherent to synaptic transmission;
Only few syts (≤4) work together for fast SV fusion on most SVs.
A high copy number of syts per SV ensures fast NT release by increasing the probability that several syts engage in Ca^{2+}/PI(4,5)P_{2} dual binding;
Binding of syts to PI(4,5)P_{2} prior to the Ca^{2+} stimulus allows some SVs to fuse very fast (submillisecond).
The molecular resolution of this model can be used to study consequences of mutations.
A sytdependent switch on the energy barrier for SV fusion
Exocytosis is a highly energydemanding reaction, for which the formation of the neuronal SNARE complex provides the energy (Jahn and Fasshauer, 2012). In our model we assume that syts regulate this process by lowering the activation energy barrier for exocytosis when they engage in Ca^{2+}/PI(4,5)P_{2} dual binding. However, how Ca^{2+}/PI(4,5)P_{2} dual binding to syt exactly reduces this energy barrier is not known. One possibility is that the energy is provided by the SNAREs themselves and that Ca^{2+}/PI(4,5)P_{2} dual binding to syt relieves a clamping function, which syt itself or the auxiliary protein complexin exerts on the SNAREs (Courtney et al., 2019; Schupp et al., 2017; Tang et al., 2006; Zhou et al., 2015; Zhou et al., 2017). Alternatively – or additionally – syt’s dual binding Ca^{2+}/PI(4,5)P_{2} might promote SNAREmediated fusion by changing the local electrostatic environment (Ruiter et al., 2019; Shao et al., 1997). Furthermore, dualbinding syts could bring SVs closer to the plasma membrane, potentially below an upper limit for full SNARE complex assembly (Araç et al., 2006; Chang et al., 2018; Honigmann et al., 2013; Hui et al., 2011; Lin et al., 2014; Nyenhuis et al., 2019; van den Bogaart et al., 2011b; Xue et al., 2008). In line with these hypothetical working mechanisms, our estimated effect each dualbound C2B domain has on the energy barrier height (~5 k_{b}T) is similar to the estimated energy barrier height for the final zippering step of the SNARE complex (Li et al., 2016). Syt’s Ca^{2+}/PI(4,5)P_{2} dual binding might also promote fusion by inducing membrane curvature or favoring lipid rearrangement (Lai et al., 2011; Martens et al., 2007). In line with this reasoning, our estimated contribution of a syt engaging in Ca^{2+}/PI(4,5)P_{2} dual binding is similar to estimates of syt1 membrane binding energies (Gruget et al., 2020; Gruget et al., 2018; Ma et al., 2017). In our model, we assume that multiple syts can simultaneously reduce the energy barrier for fusion. Here we assumed that all syts exert the same effect on this energy barrier and that the effects of more dualbound syts are additive. Whether or not this is the case will depend on the precise mechanism by which they shape the energy landscape. We show here that the simplest model (constant and independent contribution) is sufficient to reproduce the biological response.
Both the C2A and C2B domain of syt cooperate in SV exocytosis (Bowers et al., 2020; Gruget et al., 2020; Lee et al., 2013; Wu et al., 2021b). However, the exact role of the C2A domain in triggering SV fusion remains debated (FernándezChacón et al., 2002; Lee et al., 2013; Paddock et al., 2011; Sørensen et al., 2003; Stevens and Sullivan, 2003; Striegel et al., 2012). As mutation of the Ca^{2+} binding pockets of the C2A domain did not affect the affinities of Ca^{2+} and PI(4,5)P_{2} in vitro (van den Bogaart et al., 2012), we focused on the C2B domain in our model. Moreover, we aimed at developing a minimal molecular model with the least number of parameters that can fully recapitulate physical responses of the synapse. This, however, does not exclude the possibility that our C2B domainbased model indirectly describes properties of the C2A domain. For instance, Ca^{2+} binding to the C2A domain may influence the Ca^{2+} affinity of the C2B domain (Sørensen et al., 2003). This property may affect the values of other parameters of the model (e.g. our estimate of the PI(4,5)P_{2} concentration), meaning that these effects might be captured indirectly when fitting experimental data.
Allostericity buys time to synchronize syts
Our modeling study proposes that the allostericity between Ca^{2+} and PI(4,5)P_{2} binding is essential for the syts to achieve fast, synchronous, and sensitive NT release (van den Bogaart et al., 2012, Figure 4—figure supplement 1). With their experiment, van den Bogaart et al., 2012 determined steady state affinities, which do not provide information on the association/dissociation rates. This means that the allosteric effect may either be due to speeding up the association or slowing down the dissociation of Ca^{2+}/PI(4,5)P_{2} (Figure 1A). Here we implemented the latter, a reduction of the unbinding rates of both Ca^{2+} and PI(4,5)P_{2} when both species were bound to the C2B domain, which leads to a stabilization of the dualbound state. A stabilization of the Ca^{2+}bound states was also essential to reproduce the Ca^{2+} dependence of release in the previously proposed fivesite binding model (Heidelberger et al., 1994; Schneggenburger and Neher, 2000). Here we show in the context of our model that increasing the lifetime of Ca^{2+}/PI(4,5)P_{2} dual binding is particularly important to achieve fast fusion rates as it allows several C2B domains to simultaneously engage to lower the fusion barrier (Figure 4). The drawback of the strong allosteric interaction between the Ca^{2+} and PI(4,5)P_{2} bindings sites might be its potential involvement in the strong dominantnegative effects of some C2B domain mutations (Figure 6).
The stoichiometry of the SV fusion machinery
Each SV contains multiple syt copies (Takamori et al., 2006), which can jointly participate in the fusion process. However, the number of syts that can simultaneously engage with PI(4,5)P_{2} located at the plasma membrane, and thus can cooperate during fusion, is likely limited. There are several possible explanations for this limit. First, the space between the vesicular and plasma membrane is limited and crowded by many synaptic proteins (Wilhelm et al., 2014). In addition, plasma membrane association of syt may require interaction with the SNAREs (de Wit et al., 2009; Mohrmann et al., 2013; Rickman and Davletov, 2003; Zhou et al., 2015), which limits the number of association points. Moreover, the inhomogeneous distribution of PI(4,5)P_{2} in the plasma membrane might put further constraints on association of syt to PI(4,5)P_{2} (Milosevic et al., 2005; van den Bogaart et al., 2011a). Other proteins able to promote SV fusion, like Doc2, might also rely on this limited number of membrane contact points/resource and compete with syt. Our model predicts that most SVs already bind one or two slots with syt at rest (Figure 3), and this might explain the ability of syt to clamp spontaneous transmission (BouazzaArostegui et al., 2022; Courtney et al., 2019; Kochubey and Schneggenburger, 2011; Schupp et al., 2017).
We found that at least three PI(4,5)P_{2} association sites (‘slots’) were required to explain the steep Ca^{2+} dependency of neurotransmitter release (Figure 2A–C). These findings are compatible with a cryoEM analysis that identified six protein complexes between docked SVs and plasma membrane (Radhakrishnan et al., 2021). Interestingly, irrespective of the number of slots for models with three or more slots, our analysis suggests that most fusion events at [Ca^{2+}]_{i} > 1 μM occurred after engaging three syts in Ca^{2+}/PI(4,5)P_{2} dual binding (Figure 2D–E). At lower [Ca^{2+}]_{i} (0.5–1 μM, Figure 2—figure supplement 2B and Figure 2—figure supplement 3B), the number of dual bindings leading to fusion was reduced to 1–2, indicating that higher [Ca^{2+}]_{i} recruits additional syts to increase fusion rates. Although our model indicates that only few syts are involved in fusion, more syts could be involved in upstream reactions.
The predicted number of three syts involved in fast exocytosis matches experimental estimates of the number of SNAREcomplexes zippering for fast vesicle fusion (Arancillo et al., 2013; Mohrmann et al., 2010; Shi et al., 2012; Sinha et al., 2011; but higher estimates in the number of SNARE complexes actively involved in fusion have also been reported Wu et al., 2017). Moreover, our model is consistent with a previous model of neurotransmitter release at the frog neuromuscular junctions that estimated that fusion is triggered by the binding of two Ca^{2+} ions to each of three syts (Dittrich et al., 2013). That model, which describes Ca^{2+} dynamics in the AZ in detail, showed that many additional Ca^{2+} binding sites (2040) were required to enhance fusion probability, because the probability of having a single Ca^{2+} molecule in the vicinity of SVs is extremely low. Similarly, our model predicts a relevance of a high vesicular syt copy number, because, even though fusion involves only a handful of syts, many copies per SV are necessary to speed up the collision with multiple slots (Figure 5). In fact, high protein abundance could play a general role in promoting collisionlimited processes in SV fusion, and may provide an intuitive explanation for the many (~70) synaptobrevins on SVs which may assemble into SNARE complexes downstream of syt action (Takamori et al., 2006; van den Bogaart et al., 2011b).
Our model of dynamic assembly of multiple C2B domains in Ca^{2+}/PI(4,5)P_{2} dual binding in response to Ca^{2+} is fundamentally different from studies suggesting that 12–20 syts need to preassemble in higherorder complexes (rings) to execute their function in fusion (Rothman et al., 2017). A testable property to distinguish these possibilities is the sensitivity to reducing the number of syts per SV. If SV fusion relied on preassembled sytrings, it would immediately break down if the number of syts was reduced to a number preventing ring assembly, whereas our model predicts gradual effects of reduced syt copy numbers (even for titration below n_{syt} = 3; Figure 5). In line with the latter case, recent experiments that reported that SV fusion is rather sensitive to progressive reductions in syt levels (BouazzaArostegui et al., 2022). On the other hand, it might be sufficient for syts to occupy fewer slots at rest (1–2, Figure 3) to exert their effects on SV priming and clamping of spontaneous release (which are not included in our model), which might explain why these reactions appear less sensitive to syt1 reductions (BouazzaArostegui et al., 2022).
Heterogeneity in PI(4,5)P_{2} concentration between different RRP SVs
The interaction between syt and PI(4,5)P_{2} has been shown to be essential in SV exocytosis (Bai et al., 2004; Li et al., 2006; Wu et al., 2021a), but also has been found to play a role in SV docking (Chang et al., 2018; Chen et al., 2021). Consistently, we observed that at resting synapses the majority of SVs (~83%) contain at least one syt bound to PI(4,5)P_{2} in our model (Figure 3). The number of syts bound to PI(4,5)P_{2} per SV at rest highly influenced the release probability, leading to heterogeneity within the RRP (Figure 3, Wölfel et al., 2007). As PI(4,5)P_{2} levels have a large impact on release kinetics (shown in this study, but also by Walter et al., 2017), heterogeneity between RRP SVs might further be enhanced by unequal PI(4,5)P_{2} levels. Additionally, the strong impact of PI(4,5)P_{2} levels on SV fusion indicates that the dynamic regulation of PI(4,5)P_{2} occurring at the seconds time scale might strongly influence synaptic plasticity (Jensen et al., 2022).
In our model, we described PI(4,5)P_{2} levels in concentration units to constrain our model by using in vitro PI(4,5)P_{2} affinity measurements (van den Bogaart et al., 2012). However, this concentration does not only encompass the density of PI(4,5)P_{2} in the plasma membrane, but also includes the accessibility of syt to PI(4,5)P_{2}. Several studies have shown that PI(4,5)P_{2} is distributed heterogeneously over the plasma membrane in clusters that contain a high PI(4,5)P_{2} density (Honigmann et al., 2013; Milosevic et al., 2005; van den Bogaart et al., 2011a). Moreover, syts located closer to the plasma membrane will have increased access to PI(4,5)P_{2} compared to those located further away. Taken together, this indicates that the PI(4,5)P_{2} concentration is likely to vary between RRP SVs and also between individual syts on the SV. Furthermore, this implies that once a syt has engaged in PI(4,5)P_{2} binding the successive engagement of additional syts might be favored for some (those facing towards the PM) and disfavored for others (those facing from the PM). While knowledge of these details could be helpful to construct a more realistic version of our molecular model, we currently do not possess the methodology to measure these properties. Therefore, in our model, we simulated the simplest scenario where all syts have an equal probability of engaging in PI(4,5)P_{2} binding.
As the localization of syts with respect to the PM influences the accessibility of syt to PI(4,5)P_{2}, mutations in synaptic proteins and stimulation protocols that alter SV docking will affect the PI(4,5)P_{2} concentration as it is implemented in our model (Chang et al., 2018; Chen et al., 2021; Kusick et al., 2020). Using a timedependent PI(4,5)P_{2} concentration, we illustrated the impact this might have on the short term plasticity of synaptic responses (Figure 7). This is a simplification, as we did not take the individual SV/syt distances to the PM into account. This distance is affected by several synaptic proteins, including syt1, Munc13, and synaptotagmin7 (Chen et al., 2021; Imig et al., 2014; Liu et al., 2016; Quade et al., 2019; Tawfik et al., 2021; Voleti et al., 2017). A role of these proteins in shortterm plasticity is firmly established, yet precise mechanistic details are still lacking (Jackman et al., 2016; Rosenmund et al., 2002; Shin et al., 2010). The extension of models based on molecular interactions such as presented here should allow reproduction of responses to more complex synaptic activity patterns relevant for neural processing. Particularly the molecular resolution of such models will be useful to conceptualize the importance of specific molecular interactions for physiological and pathological processes at the synapse.
Materials and methods
In this paper, we propose a model for SV fusion induced by Ca^{2+} and PI(4,5)P_{2} binding to n_{syts} syts per SV, with at most M_{slots} syts per SV engaging in PI(4,5)P_{2} binding at the same time. We implemented the model in two ways for different simulation purposes: (1) an implementation based on the analytical solution of the model, and (2) an implementation following the Gillespie algorithm (Gillespie, 2007) (Matlab procedures for simulations can be found in Source code 1). In the first implementation, we assume a constant [Ca^{2+}]_{i,} (allowing us to simulate Ca^{2+} uncaging experiments: SV reactions at a Ca^{2+} level reached by the uncaging stimulus from a steady state starting point calculated for the resting Ca^{2+} concentration of 50 nM), whereas the second version was implemented to allow for [Ca^{2+}]_{i} to vary over time (allowing us to simulate APevoked responses). Another important difference between the two approaches is that the analytical solution describes the binding state of an entire SV and the Gillespie version describes the binding state of each individual syt. Both implementations allow for stochastic evaluation of the model. The first implementation is used in Figure 2, Figure 2—figure supplement 4, Figure 2—figure supplement 5, Figure 2—figure supplement 6, Figure 2—figure supplement 7, Figure 3, Figure 3—figure supplement 1, Figure 4, Figure 4—figure supplement 1, Figure 5, Figure 5—figure supplement 1, Figure 5—figure supplement 2. The second implementation is used to simulate the APevoked responses and individual SV binding states in Figure 2—figure supplement 2, Figure 2—figure supplement 3, Figure 2—figure supplement 4, Figure 2—figure supplement 5, Figure 2—figure supplement 6 Figure 4, Figure 5, Figure 5—figure supplement 1, Figure 5—figure supplement 2 Figure 6, and Figure 6—figure supplement 1, and Figure 7. Consistency between the two approaches was validated by comparison of simulation result distributions in quantilequantile (QQ) plots (Figure 5—figure supplement 1).
SV states and possible reactions in the analytical version of the model
Request a detailed protocolIn the analytical solution of the model, we describe for each SV the number of syts having bound two Ca^{2+} ions, PI(4,5)P_{2}, or both species. Since syts were assumed to work independently, their order is not relevant, and we therefore do not need to describe the binding state of each individual syt. The possible binding states of an SV are described in Figure 1—figure supplement 2. Each state is represented by the triplet (n,m,k), with k denoting the number syts having bound PI(4,5)P_{2}, m denoting the number of syts having bound two Ca^{2+} ions, and n denoting the number of syts having bound both species and thereby having formed a dual binding. n_{syts} is the total number of syts per SV. M_{slots} restricts the number of syts having bound PI(4,5)P_{2} simultaneously, which includes syts having bound PI(4,5)P_{2} only (k) and those having formed a dual binding (n). Taken together, this implies that for all states in the model, it holds that
We numbered the states systematically following a lexicographic ordering, excluding the states that violate the inequalities described above. To illustrate, we write the ordering of all the states (m,n,k) with n_{syts} = 3 and M_{slots} = 2:
(0,0,0),(0,0,1),(0,0,2),(0,1,0),(0,1,1),(0,1,2),(0,2,0),(0,2,1),(0,3,0),(1,0,0),(1,0,1),(1,1,0),(1,1,1),(1,2,0),(2,0,0),(2,1,0).
Besides these binding states, an additional state, F, describes whether the SV has fused. With n_{syts} = 15 and M_{slots} = 3, a single SV in our model has 140+1 states. From any state, there are at most 9 possible reactions, one being SV fusion and the other 8 being (un)binding of Ca^{2+} or PI(4,5)P_{2} to/from a syt. The rates for the possible reactions of a single SV in this model are summarized in Table 3. In many cases, only a subset of the 8 (un)binding reactions are allowed because of the inequalities above (noted under ‘Condition’ in the table).
The reaction rates of (un)binding Ca^{2+} or PI(4,5)P_{2} are calculated as the number of syts available for (un)binding (computed using n_{syts}, n, m, k) times the reaction rate constant (α, β, γ, δ), and, in the case of binding reactions, times the concentration of the ligand ([PI(4,5)P_{2}] or [Ca^{2+}]). We assumed binding of two Ca^{2+} ions to a single C2B domain. In our model, this twostep process is simplified to a single reaction step by taking [Ca^{2+}]_{i} to the power of two. This simplification is reasonable, because we assumed that syt could only associate to the vesicular membrane when two Ca^{2+} ions are bound, and binding of one Ca^{2+} ion would not induce an ‘intermediate’ association state to the membrane, nor would it affect the allosteric interaction. By simplifying the two Ca^{2+} binding/ubinding steps to one, we indirectly assumed high cooperativity between the two Ca^{2+} binding sites. To account for the limit on the number of syts bound to PI(4,5)P_{2}, the number of available, empty slots, (M_{slots}nk), was multiplied on the PI(4,5)P_{2} binding rates. The fusion rate of the SV is computed by L_{+}·f^{n} (similar to Lou et al., 2005), with L_{+}denoting the basal fusion rate and f the factor of increase in fusion rate by each dual binding being formed. L_{+}was set to 4.23e4 s^{–1} to match the release rate measured at low [Ca^{2+}]_{i}, given an average size of the RRP of 4000 SVs (see below).
The affinities for Ca^{2+} and PI(4,5)P_{2} binding to syt were set to previously determined dissociation constants (K_{D,Ca2+} = β/α = 221^{2} µM^{2}, K_{D,PI(4,5)P2} = δ/γ = 20 µM) obtained using in vitro microscale thermophoresis experiments (van den Bogaart et al., 2012). For determination of the dissociation constant of Ca^{2+}, van den Bogaart and colleagues assumed binding of a single Ca^{2+} ion to the C2AB domain (van den Bogaart et al., 2012). The K_{D,2Ca2+} of the reaction describing binding of two Ca^{2+} ions, in our model, can be computed from the experimentally derived dissociation constant by taking it to the power of two. This was corroborated by refitting the experimental data with a hill coefficient of 2, which yielded a similar K_{D,2Ca2+} value of~221^{2} (data not shown).
The in vitro experiments revealed a change in syt1 Ca^{2+} affinity upon binding PI(4,5)P_{2}, and vice versa (van den Bogaart et al., 2012), indicating a positive allosteric relationship between the two species. We assumed this allosteric effect was due to a stabilization of the dualbound state by lowering of the unbinding rates of Ca^{2+} and PI(4,5)P_{2} with a factor (A=(3.3/221)^{2} = 0.00022) and occurs when both species have bound. Upon dual binding, both rate constants for unbinding Ca^{2+} and PI(4,5)P_{2} are multiplied by A, since any closed chemical system must obey microscopic reversibility (Colquhoun et al., 2004). Using the biochemically defined affinities, the number of free parameters in our model was constrained to:
The values of β and δ were determined according to the affinities for each choice of α and γ.
The steady state of the system
Request a detailed protocolThe steady state of the system before stimulation was determined at a resting, global [Ca^{2+}]_{i} of 0.05 µM (except for simulations with Ca^{2+} levels below this basal value, for those we assumed [Ca^{2+}]_{rest}=[Ca^{2+}]_{i}). To compute the steady state, we assumed that no fusion took place, ignoring the very low fusion rate at resting [Ca^{2+}]_{i}. Under these conditions, the model is a closed system of recurrent states and obeys microscopic reversibility, that is for every closed loop state diagram, the product of the rate constants around the loop is the same in both directions (Colquhoun et al., 2004). Microscopic reversibility implies detailed balance, meaning that every reaction is in equilibrium at steady state. Thus, for any two states S_{i} and S_{j} which are connected by a reaction, the steady state distribution obeys
where [S_{i}] and [S_{j}] are steady state quantities and r_{ij} and r_{ji} are the reaction rates between S_{i} and S_{j}. Using this property, we calculated the steady state iteratively by setting the population of the first state (state (0,0,0)) to 1, and thereafter iteratively computing the population of the following state (following the lexicographic ordering as described above) using the following formulae:
Afterwards, each state was divided by the sum of all state values and multiplied by the number of SVs in the RRP. In our model simulations, the size of the RRP was variable and followed a gamma distribution with a mean of 4000 SVs and a standard deviation of 2000 SVs (see Figure 2—figure supplement 1), based on experimental estimates from the calyx of Held (Wölfel and Schneggenburger, 2003). In the following calculations, we use $\phi$ to denote the steady state probability vector (i.e. normalised to sum to 1).
Computation of fusion probabilities and fusion rate
Request a detailed protocolThe analytical implementation of our model allowed us to compute the fusion rate and cumulative fusion probabilities with a constant [Ca^{2+}]_{i} after stimulus onset (t=0), thereby mimicking conditions in Ca^{2+} uncaging experiments. The constant [Ca^{2+}]_{i} makes the model a homogenous Markov Model. The transition rates of the model can be organized in the intensity matrix, Q, such that,
where k_{i,j} is the rate of the reaction from state i to state j. Given initial conditions, φ (steady state normalized to a probability vector) and intensity matrix Q_{ξ,C} ( ξ being the free model parameters and C being the Ca^{2+} concentration),
is the distribution of SV states at time t, i.e. a 1 x n_{states} vector with element i being the probability of being in state i at time point t. The single SV cumulative fusion probability (being a function, G, of time, t, defined by ξ and C) is the last element, which we will denote with a subscript F,
The fusion rate of a single SV can be calculated directly as the last element of the derivative of (Equation 1):
Multiplying (Equation 1) and (Equation 2) with the number of SVs yields the cumulative fusion function and fusion rate function, respectively. For simulation of release rates and release latencies (in Figures 2—5) and some figure supplements, we computed (Equation 1) and (Equation 2) using the best fit parameters from fitting with M_{slots} = 3 and n_{syts} = 15. This was done for 31 [Ca^{2+}]_{i} values ranging from 0.001 µM to 80 µM (0.001, 0.1, 0.2, …., 0.9, 1, 1.25, 1.5, 1.75, 2, 2.5, 3, 4, …, 9, 10, 20, 30, …, 80 µM) for t∈[0,100] ms with a time step of 0.01ms. In addition, the functions were calculated in the same way using the best fit parameters for M_{slots} = 1,2,4,5,6 with n_{syts} = 15 for simulations depicted in (Figure 2) and figure supplemets. In some conditions, especially at low [Ca^{2+}]_{i}, a longer span of the cumulative fusion probability was required and was calculated with the same time step size.
Computation of peak release rates
Request a detailed protocolThe peak of the fusion rate can be computed by multiplying the maximum value of the single SV fusion rate function, (Equation 2), with n_{ves}. To allow for a variable RRP size, a set of 1000 n_{ves} values were drawn according to the RRP size distribution, the peak release rates were determined, and the mean and 95% prediction interval determined for each Ca^{2+} concentration.
For parameter exploration (Figure 2G) and for computing the release rates in the fitting routine, it was not feasible to calculate the fusion rate over 100ms with high temporal precision as described above. Instead, we implemented a custom search algorithm (scripts can be found in accompanying zipfile “Source_code1.zip”), which was constructed to shorten calculation time by taking advantage of the release rate function being unimodal. We first found a time point, t_{max}, at which 75–90% of the SVs had fused. Having computed different time points in the time interval [0,t_{max}], gave us an interval in which the fusion rate showed a local maximum. The algorithm then narrowed the time interval down until a time of peak was found with a precision of 0.01ms. This method shortened simulation time considerably.
Stochastic simulation of release latencies
Request a detailed protocolRelease latencies, which are defined as the time point of the fifth SV fusion event after the onset of simulation, were simulated stochastically by drawing p_{i}∈(0,1), i=1,…,n_{ves}, from the uniform distribution for each of the 1000 repetitions and each evaluated [Ca^{2+}]_{i}. Each of these random numbers corresponds to the fusion time of an SV, which can be determined interpolation of the single SV cumulative fusion probability function (Equation 1). To obtain the time point of the fifth SV fusion, the fifth lowest p_{i} was used. In the corresponding figures, the medians were plotted, since the probability distribution of the release latencies (derived below) was skewed, and the reported data points were single measurements.
Fitting the model to experimental data
Request a detailed protocolWe next fitted the model to already published data describing the Ca^{2+} dependence of NT release in the mouse calyx of Held (Kochubey and Schneggenburger, 2011). The data consist of measurements from Ca^{2+} uncaging experiments, where the release latency, defined as the time point of the fifth SV fusion event, and the peak release rate were estimated at different [Ca^{2+}]_{i}. Besides the four free model parameters, $\xi $, an additional parameter, d, was fitted. d is a constant added to the release latencies computed from the model to account for the experimentally observed delay (Kochubey and Schneggenburger, 2011).
Since the variance in the experimental data points also contains information on the underlying biological mechanism, we wanted to take the distribution of individual data points into account when obtaining estimates of the unknown parameters. We therefore derived the likelihood function, which describes how well the model captures the distribution of the release latencies. Obtaining this function for the peak release rates was not feasible. The experimental peak release rates were therefore compared to the average model prediction. Both measures of describing the correspondence between model simulations and experimental data were combined in a cost value which was optimized to estimate the best fit parameters (the lower this cost value the better the correspondence between model predictions and experimental data).
The best fit was obtained by minimizing the following cost function:
where i=1,…, ${n}_{rates}$ and j=1,…, ${n}_{latencies}$ are the indices of the experimental [Ca^{2+}]_{i},
are the squared deviations of the peak release rates (1 /ms^{2}) and $\ell$ is the logarithm of the likelihood of the release latencies (see derivation below). To combine the two measures of distance between model and experimental data, the squared deviation of the peak release rates was multiplied by a factor 2 before subtracting the logarithm of the likelihood of the release latencies. The cost value was minimized using the inbuilt Matlab function fminsearch, which uses a NelderMead simplex method. fminsearch was run with different initial parameter values to verify that the global minimum of the cost function was found. During the fitting, the lower and upper bounds of d were set to, respectively, 0.3ms and 0.405ms (with the upper bound corresponding to the smallest release latency in the experimental data set). α,γ, and [PI(4,5)P_{2}] had an upper bound of 10^{10}, and all free parameters needed to be positive. The maximum number of iterations and function evaluations was set to 5000.
The likelihood function of release latencies with fixed RRP size
Request a detailed protocolTo fit the model to the experimental release latency measurements, we derived the likelihood function, which is the probability density function of the model for given parameters evaluated at the experimental data points. Thus, optimizing the likelihood function yields parameters for which the data points are most likely if the model is true. We first derive the likelihood of release latency in the case of a fixed RRP size (n_{ves}).
We define the stochastic variable X, which describes the stochastic process of the state of a single SV in our model. The fusion time of the SV, τ, is defined as
where τ itself is a stochastic variable. We define a stochastic vector, Z, which consists of all n_{ves} fusion time points in a single experiment. They come from independent, identically distributed stochastic processes with cumulative distribution function G_{ξ,C}(t), given in (Equation 1). As the release latency is defined as the time of the fifth SV fusion, we order the Z variable outcomes (z_{(1)}, z_{(2)}, …, z_{(nves)}) from first to last fusion time. Using the transformation
we obtain a sequence of stochastic variables, U_{(i)}, which are uniformly distributed on the interval (0,1). The ordering is preserved, since G_{ξ,C}(t) is monotonically increasing, and U has probability density function
with respect to t. G_{ξ,C}’(t)≥0 is given in (Equation 2). From order statistics it follows that the k^{th} ordered U is betadistributed with probability density
where u=G_{ξ,C}(t). Thus, the transformed variable U_{(k)}, is betadistributed, with
In the case of the release latency, we are interested in the fifth fusion event (k=5). Thus, with a fixed RRP size, the likelihood value for the release latency observations, T^{*}=(t^{*}_{1}, t^{*}_{2}, …, t^{*}_{M}), at all M Ca^{2+} concentrations is
In the optimization we minimize minus the loglikelihood:
which is equivalent to maximizing the likelihood function.
The likelihood of release latencies with variable RRP size
Request a detailed protocolIn our model, the RRP size is assumed to follow a Gamma distribution. Let x denote the RRP size, f_{RRP}(x) the probability density of the Gamma distribution, and u=G_{ξ,C}(t) as defined above. The probability density of the release latency at Ca^{2+} concentration C_{i} is given by
where K is a normalization constant, K=1P(x<5)≈1. The lower limit of the integral reflects that the release latency is only defined when there are more than 5 SVs in the RRP. In simulations this corresponds to redrawing the RRP size whenever an RRP size <5 SVs occurs, which happens with probability ~3e11, and is accounted for in the normalization constant K in the following. Inserting the probability density function of a Gamma distribution with shape parameter k and scale parameter θ, we get:
We now define the following variables:
By factoring out and substituting in the above equation we get
Furthermore, we have
Since
we can derive the following useful formula:
The third equality follows from the fact that the function in the second integral from above is the probability density function of a gamma distribution with shape parameter k+n and scale parameter $\stackrel{~}{\theta}$ . We therefore replace it with the cumulative distribution function of the gamma distribution, where $\gamma $ is the lower incomplete gamma function. Combining (Equations 4–7) yields an explicit expression of the likelihood of a single delay in (Equation 3), as
with
We then minimize the sum of minus the loglikelihoods of the release latency observations.
Syt binding states in the Gillespie simulation of model
Request a detailed protocolIn the Gillespie algorithm, the binding state of each individual syt is tracked. The state of the system at time point t, X(t), is given by a n_{syt} x n_{ves} matrix. Each element in this matrix describes the binding state of an individual syt using a twodigit coding system; 00 for no species bound to syt, 01 for PI(4,5)P_{2} bound, 10 for two Ca^{2+} ions bound, and 11 for both species bound (dualbinding syt). As with the analytical implementation, each syt can undergo a subset of the 8 different (un)binding reactions (Figure 1A), depending on the binding state of the respective syt. The fusion rate, which depends on the number of dualbound syts per SV, is determined for the entire SV.
Determining the initial state of the system
Request a detailed protocolThe steady state (initial state, X(0)) was computed using the same method as described above (see section ‘The steady state of the system’) using [Ca^{2+}]_{i} = 0.05 µM as the resting condition. This resulted in φ, the probability vector of a single SV to be in the different SV states at steady state. To stochastically determine X(0), we first determined the binding state for each SV, that is how many dual bindings are formed (n) and how many syts have bound Ca^{2+} (m) and how many PI(4,5)P_{2} (k). For that we drew p_{j}∈(0,1), j=1…,n_{ves}, from the uniform distribution. The state number of the j^{th} SV, s, was determined by:
Via the ordering of states explained above, s can be linked to the state triplet (n_{s},m_{s},k_{s}). As the order of syts is irrelevant for model simulation, this information on the state of SV_{j} was transferred to the j^{th} column of the X(0) matrix in a systematic way: The first n_{s} elements were labeled with ‘11’; elements n_{s} +1 to n_{s} +m_{s} were labeled with ‘10’; and elements n_{s} +m_{s} + 1 to n_{s} +m_{s} + k_{s} were labeled with ‘01’, and the remaining elements (n_{s} +m_{s} + k_{s} +1) to n_{syt} were set to ‘00’.
Gillespie algorithmbased simulations of the model
Request a detailed protocolFor stochastic evaluation of the model by the Gillespie algorithm (Gillespie, 2007), we next introduced the propensity function B, which is defined by:
${B}_{i,j}\left(x\right)dt\text{}:=$ the probability given (X(t)=X), that the i^{th} syt of the j^{th} SV will undergo a reaction in the next infinitesimal time interval [t, t+dt].
For element i,j in B, the total reaction propensity is the sum of propensities of possible reactions given the binding state of the syt and can be computed as follow:
with PIP_{tot,j} the total number of syts of SV j bound to PI(4,5)P_{2}. To include the propensity for fusion of an SV, an additional row in B, which we index with the denotation f, describes the propensity of fusion for each SV in the matrix;
This makes B a matrix of (n_{syt} +1)×n_{ves}. We denote the sum of all elements in B by B_{0}. Using B_{0} and 3 random numbers (r_{n}∈(0,1), n=1,2,3) drawn from the uniform distribution, we determined the time step to the next reaction and which SV and syt this reaction affects. The time to next reaction, $\stackrel{~}{\tau}$, is given by
since it is exponentially distributed with rate B_{0}. The index, j, of the SV undergoing a reaction is the first index that satisfies:
Similarly, the index i of the syt in SV j undergoing a reaction is the smallest integer fulfilling:
If the row index i equals n_{syt} +1, a fusion reaction occurs. The fusion time (t + $\stackrel{~}{\tau}$) is saved and all the propensities of SV j in B are set to 0. If i is smaller or equal to n_{syt} a binding or unbinding reaction of one of the two species occurs. To determine which of the four possible reactions is occurring, we define an additional propensity vector, b_{react}. The first element in b_{react} denotes the propensity of PI(4,5)P_{2} binding, the second element the propensity of Ca^{2+} binding, and the third and fourth element the unbinding of PI(4,5)P_{2} and Ca^{2+}, respectively. These elements are given by:
Additionally, we define the transition matrix V, which describes the change in the state of X_{i,j} induced by the four reactions:
A fourth random number, r_{4}∈(0,1), drawn from the uniform distribution, determines which reaction, h, occurs:
The state of the corresponding SV and syt, X_{i,j}, is replaced by ${X}_{i,j}+V{}_{h}$ and t by $t+\tau $. Then B is updated according to the change in X, and all steps are repeated. This iterative process continues until all SVs are fused, when simulating the model with a fixed [Ca^{2+}]_{i}.
When simulating APevoked responses (Figures 5 and 6), we used a Ca^{2+} transient describing the microdomain [Ca^{2+}]_{i} sensed locally by primed SVs in the mouse calyx of Held upon AP stimulation (Wang et al., 2008). This Ca^{2+} transient also formed the basis for the Ca^{2+} signal used to simulate a paired pulse stimulus (Figure 7), where the transients were placed with a 10 ms interval. Additionally, for the paired pulse stimulus, we added a residual Ca^{2+} transient to the signal (exponential decay with amplitude: 0.4 µM, decay time constant: 0.154 s^{–1}). Similar to the uncaging simulations, the [Ca^{2+}]_{i} before the onset of the stimulus was 0.05 µM. Since the Ca^{2+} concentration is a factor in the reaction rates, the propensity matrices B and b_{react} were not only updated to the new state matrix, $X\left(t+\stackrel{~}{\tau}\right)$, but also to a new [Ca^{2+}]_{i} at time $t+\stackrel{~}{\tau}$. [Ca^{2+}]_{i} at time $t+\stackrel{~}{\tau}$ was determined from the Ca^{2+} transient using linear interpolation, and B and b_{react} were updated correspondingly. As the propensity of Ca^{2+} binding is largely dependent on [Ca^{2+}]_{i}, the time step between two updates in [Ca^{2+}] and the propensity matrices was set to be at most 8e^{–4}s. If $\stackrel{~}{\tau}$ determined from B at time t was larger than this time step, no reaction occurred, but the system and [Ca^{2+}] were updated to time t+0.8ms. The model was evaluated until the end of the Ca^{2+} transient. Similarly, in the case of a variable PI(4,5)P_{2} transient (Figure 7), the PI(4,5)P_{2} concentration was updated at least every 0.8 ms. Because this approach requires simulation of all individual (un)binding reactions and fusion events it is not feasible to perform 1000 repetitions. Instead, simulations were repeated 200 times (Figures 5 and 6). Like with the computation of the release latencies and maximal fusion rates, we assumed a variable RRP. However, instead of drawing random pool sizes from a gamma distribution, we used the 200 quantiles of the pdf of the RRP sizes, because of the limited number of repetitions and the large impact of the RRP size on the model predictions. In Figure 7, we reduced the number of repetitions to 50, to represent an experimental condition more closely. The set of RRP sizes was drawn randomly once. This set of random RRP values was used for both the mutant and the WT condition displayed in the figure.
Simulating the model with mutant syts
Request a detailed protocolFor mutations in syt that affect the binding and unbinding rates of PI(4,5)P_{2} and Ca^{2+}, the procedure described above was repeated with adjusted parameters when simulating a model containing only mutant syts. For a model in which mutant proteins were expressed together with WT syts (simulations of heterozygous condition in Figure 6), the procedure was changed slightly.
For a model with p WT and q mutant syts, the number of states of an SV increases drastically and was now described by six values; the number of WT syts bound to Ca^{2+}, PI(4,5)P_{2} or both, and the number of mutant syts bound to Ca^{2+}, PI(4,5)P_{2} or both. The dimensions of the Qmatrix used to compute the steadystate probability of a single SV increased with n_{states}. X(t=0) was computed using the same principle as described above, with the important difference that the first p rows represented the binding status of the WT syts, and row p+1 to p+q that of the mutants. In B this ordering of WT and mutant syts is the same. The parameters used to compute b_{react} depended on whether a reaction occurred to a WT syt, $i\le p$, or a mutant syt, ${n}_{syt}\ge i>p$.
Simulation of EPSCs
Request a detailed protocolSimulated EPSCs were obtained using both model implementations. The analytical implementation of our model was used to simulate fusion times for a constant [Ca^{2+}]_{i} (Figure 4D and G). The Gillespie version of the model was used to simulate APevoked EPSCs with or without mutant syts (Figure 5C–E, Figure 5—figure supplement 2EF, Figure 6C–D, Figure 6—figure supplement 1, and Figure 7B–D). In both approaches, the stochastically determined fusion times, determined as described above, were rounded up to the next 0.02ms, leading to a sampling rate of 50 kHz. The sampled fusion times were convolved with a mEPSC to generate simulated EPSCs. The standard mEPSC used for deconvolution followed the equation described by Neher and Sakaba, 2001:
with τ_{1}=0.12 ms (time constant of fast decay), τ_{2}=13 ms (time constant of slow decay), τ_{0}=0.12 ms (time constant of rise phase), ρ=1.e5 (proportion of slow phase in decay), and A being a normalization constant making the amplitude 60 pA. Parameter values were chosen such that the kinetics of the mEPSC would match events measured in the Calyx of Held (Chang et al., 2015). In Figure 4D and G traces show three randomly chosen eEPSCs in each panel. Representative eEPSC traces shown in Figure 5, Figure 5—figure supplement 2, Figure 6, Figure 6—figure supplement 1, and Figure 7 are simulated eEPSCs with the amplitude closest to the mean eEPSC amplitude.
Simulating APevoked EPSCs with variable number of syt
Request a detailed protocolTo investigate the effect of variability in the number of syts expressed per SV on variance between simulated APevoked traces (Figure 5), we first had to determine the steady state. For this we computed the probability vector of a single SV to be in the different SVstates at steady state (φ) for n_{syt} = 1,…,50. Subsequently, for each SV and each repetition (n=1000) a random number of n_{syt} was drawn from the Poisson distribution with mean = 15. When the value 0 was drawn, it was replaced by n_{syt} = 1. Using these values and φ determined for each n_{syt}, we computed the steady state matrix (X(0)) as described above (‘Determining the initial state of the system’). To reduce computation time, we evaluated a model containing 100 vesicles 40 times instead of evaluating 4000 vesicles simultaneously. The fusion times obtained when driving the model with the Ca^{2+} transient were combined afterwards. This is valid since all SVs act independently in the model. For the condition with a variable RRP size, fusion times were selected randomly until the RRP size of that specific repetition was reached. Afterwards, the fusion times were convolved with the mEPSC to obtain simulated APevoked responses. To quantify the variance in the traces (Figure 5E), we computed the standard deviation of the simulated eEPSCs at each data point (300 data points corresponding to the time interval 0–6ms) and summed those values.
Data availability
All data and software codes generated and used during this study are included in the manuscript and supporting files. Source data is included for all figures.
References

Close membranemembrane proximity induced by ca(2+)dependent multivalent binding of synaptotagmin1 to phospholipidsNature Structural & Molecular Biology 13:209–217.https://doi.org/10.1038/nsmb1056

Titration of syntaxin1 in mammalian synapses reveals multiple roles in vesicle docking, priming, and release probabilityThe Journal of Neuroscience 33:16698–16714.https://doi.org/10.1523/JNEUROSCI.018713.2013

PIP2 increases the speed of response of synaptotagmin and steers its membranepenetration activity toward the plasma membraneNature Structural & Molecular Biology 11:36–44.https://doi.org/10.1038/nsmb709

Ca(2+)synaptotagmin directly regulates tSNARE function during reconstituted membrane fusionNature Structural & Molecular Biology 13:323–330.https://doi.org/10.1038/nsmb1076

Deconstructing synaptotagmin1’s distinct roles in synaptic vesicle priming and neurotransmitter releaseThe Journal of Neuroscience 42:2856–2871.https://doi.org/10.1523/JNEUROSCI.194521.2022

Dynamic binding mode of a synaptotagmin1SNARE complex in solutionNature Structural & Molecular Biology 22:555–564.https://doi.org/10.1038/nsmb.3035

How to impose microscopic reversibility in complex reaction mechanismsBiophysical Journal 86:3510–3518.https://doi.org/10.1529/biophysj.103.038679

Structure/function analysis of ca2+ binding to the C2A domain of synaptotagmin 1The Journal of Neuroscience 22:8438–8446.https://doi.org/10.1523/JNEUROSCI.221908438.2002

Stochastic simulation of chemical kineticsAnnual Review of Physical Chemistry 58:35–55.https://doi.org/10.1146/annurev.physchem.58.032806.104637

Phosphatidylinositol 4,5bisphosphate clusters act as molecular beacons for vesicle recruitmentNature Structural & Molecular Biology 20:679–686.https://doi.org/10.1038/nsmb.2570

Mechanism and function of synaptotagminmediated membrane appositionNature Structural & Molecular Biology 18:813–821.https://doi.org/10.1038/nsmb.2075

Biophysical physiology of phosphoinositide rapid dynamics and regulation in living cellsThe Journal of General Physiology 154:e202113074.https://doi.org/10.1085/jgp.202113074

The readily releasable pool of synaptic vesiclesCurrent Opinion in Neurobiology 43:63–70.https://doi.org/10.1016/j.conb.2016.12.012

A postdocking role of synaptotagmin 1C2B domain bottom residues R398/399 in mouse chromaffin cellsThe Journal of Neuroscience 35:14172–14182.https://doi.org/10.1523/JNEUROSCI.191115.2015

Synaptic vesicles transiently dock to refill release sitesNature Neuroscience 23:1329–1338.https://doi.org/10.1038/s41593020007161

Synaptotagmin 1 modulates lipid acyl chain order in lipid bilayers by demixing phosphatidylserineThe Journal of Biological Chemistry 286:25291–25300.https://doi.org/10.1074/jbc.M111.258848

Phosphatidylinositol phosphates as coactivators of ca2+ binding to C2 domains of synaptotagmin 1The Journal of Biological Chemistry 281:15845–15852.https://doi.org/10.1074/jbc.M600888200

Mutations in the second C2 domain of synaptotagmin disrupt synaptic transmission at Drosophila neuromuscular junctionsThe Journal of Comparative Neurology 436:4–16.

How synaptotagmin promotes membrane fusionScience 316:1205–1208.https://doi.org/10.1126/science.1142614

Synaptotagmin interaction with SNAP25 governs vesicle docking, priming, and fusion triggeringThe Journal of Neuroscience 33:14417–14430.https://doi.org/10.1523/JNEUROSCI.123613.2013

Combining deconvolution and noise analysis for the estimation of transmitter release rates at the calyx of heldThe Journal of Neuroscience 21:444–461.

Dual roles of the C2B domain of synaptotagmin I in synchronizing ca2+dependent neurotransmitter releaseThe Journal of Neuroscience 24:8542–8550.https://doi.org/10.1523/JNEUROSCI.254504.2004

Membrane penetration by synaptotagmin is required for coupling calcium binding to vesicle fusion in vivoThe Journal of Neuroscience 31:2248–2257.https://doi.org/10.1523/JNEUROSCI.315309.2011

The ca2+ affinity of synaptotagmin 1 is markedly increased by a specific interaction of its C2B domain with phosphatidylinositol 4,5bisphosphateThe Journal of Biological Chemistry 284:25749–25760.https://doi.org/10.1074/jbc.M109.042499

Mechanism of calciumindependent synaptotagmin binding to target snaresThe Journal of Biological Chemistry 278:5501–5504.https://doi.org/10.1074/jbc.C200692200

Munc13 C2B domain is an activitydependent ca2+ regulator of synaptic exocytosisNature Structural & Molecular Biology 17:280–288.https://doi.org/10.1038/nsmb.1758

Calciumdependent docking of synaptic vesiclesTrends in Neurosciences 44:579–592.https://doi.org/10.1016/j.tins.2021.04.003

Examining synaptotagmin 1 function in dense core vesicle exocytosis under direct control of ca2+The Journal of General Physiology 122:265–276.https://doi.org/10.1085/jgp.200308855

Synaptotagmins: why so many?The Journal of Biological Chemistry 277:7629–7632.https://doi.org/10.1074/jbc.R100052200

Synaptotagmin1 may be a distance regulator acting upstream of SNARE nucleationNature Structural & Molecular Biology 18:805–812.https://doi.org/10.1038/nsmb.2061

Phosphatidylinositol 4,5bisphosphate increases ca2+ affinity of synaptotagmin1 by 40foldThe Journal of Biological Chemistry 287:16447–16453.https://doi.org/10.1074/jbc.M112.343418

The janusfaced nature of the c(2)b domain is fundamental for synaptotagmin1 functionNature Structural & Molecular Biology 15:1160–1168.https://doi.org/10.1038/nsmb.1508
Decision letter

Henrique von GersdorffReviewing Editor; Oregon Health and Science University, United States

Richard W AldrichSenior Editor; The University of Texas at Austin, United States

Victor MatveevReviewer; New Jersey Institute of Technology, United States
Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.
Decision letter after peer review:
Thank you for submitting your article "Allosteric stabilization of calcium and lipid binding engages three synaptotagmins in fast exocytosis" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Richard Aldrich as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Victor Matveev (Reviewer #1).
The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.
Essential revisions:
1) Lines 670672: The main assumption of the model is that nearsimultaneous binding of 2 calcium ions greatly increases the affinity of PIP2 binding. However, simultaneous binding of 2 Ca^{2+} ions must imply an even stronger cooperativity between the two Ca^{2+} binding events than the cooperativity between Ca^{2+} and PIP2 binding. While equilibrium properties are not affected by such an assumption, the timedependence of the model should depend on the precise sequence of binding events. Would the release latency predictions change if both Ca^{2+} binding events were resolved? Of course, the stochastic simulations for a more complicated model would be costly, and extra model parameters are not desired, but one trial analytic solution for a constant Ca^{2+} level is probably not hard. I am not insisting that any new simulations be included in the paper, but it would help to discuss briefly whether the assumption of two simultaneous Ca^{2+} binding affects the release latency estimation.
2) On page 10, the Authors describe their checks to ensure that the model parameter optimization procedure is not stuck at a local minimum. However, a more straightforward common check would involve repeating the fminsearch parameter optimization algorithm multiple times, while starting at different initial values of parameters. Have the Authors performed such a check? This would essentially quantify uncertainties in inferred parameter values in Table 2.
3) In the discussion (Line 542) the authors note that their model is fundamentally different from the proposed preassembly of Syt rings at the base of a vesicle, and that this is big testable prediction. How could this be tested? And do Syts actually preassemble with slots in their model? Figure 3B shows the number of crosslinks that precede fusion, and it would be nice to see how many Syts are bound to PIP2 at rest for any given number of "slots".
4) I like Figure 2 Sup2, where the authors probe the number of crosslinks prior to fusion. Interestingly, the speed of Ca ramping affects the number of crosslinks, and a significant number of events occur after 4 crosslinks only when Ca increases at the fastest rate (0.001 s). This is still much slower than the expected rate of Ca nanodomain formation (<100 us). Why did the authors explore only extremely slow Ca ramp rates? It seems that the ramping rates of 10 us to 10 ms would more accurately affect the condition of action potential induced Ca influx.
5) It is not immediately clear to me why the model with 6 slots never predicts fusion with 56 crosslinks. Why does the model fail to fit well if each crosslink contributes less to lowering the energy barrier, and thus 56 crosslinks are required to drive fusion?
6) If a single crosslink brings the vesicle closer to the PM, this will introduce a new form of allostery by increasing the effective concentration of PIP2 sensed by Syts. Is this reflected in the model?
7) In the model, PIP2 binding is required for Syt1driven fusion. As a result, only Syts with a "slot" can participate in fusion. However, while K325, 327A mutations that disrupt Caindependent PIP2 binding lead to less synchronous release, synchronous release can be restored using a pairedpulse stimulation paradigm where an initial action potential drives Cadependent membrane attachment of vesicles to the plasma membrane. This suggests that PIP2 binding is not required for "cross linking". How might the model change if Syts without slots could also contribute to fusion in your model?
8) The evaluation of mutants in Figure 5 should be better tied to actual biology. What mutation does "Ca^{2+} binding", and "Aon" mimic? No citations are provided for the many studies where Syt mutants were expressed in Syt1 KO neurons. This seems like a perfect place to test the role of PIP2 binding with more complicated Ca stimuli. Could the model be adapted to explore the effect of K/A mutations described by Chang, Trimbuch and Rosenmund?
9) Recent experimental and modeling work with Synaptotagmin and SNAREs has been published (Wu et al., 2021; eLife; https://elifesciences.org/articles/68215). These authors say: "To test whether Syt1 affected fusion pores in this system, we coreconstituted ~4 copies of recombinant fulllength Syt1 together with ~4 copies of VAMP2 (per disc face) into large nanodiscs called nanolipoprotein particles…". The work seems particularly relevant to your paper. Also please take a look at another paper from this group that suggests fast exocytosis requires up to 15 SNAREs and Syt complexes (Wu et al., eLife, 2017; https://elifesciences.org/articles/22964). Please read and discuss in your paper these recent modeling studies, which seem to suggest that large numbers of Syt's and SNAREs are needed for fast exocytosis at synapses.
Reviewer #1 (Recommendations for the authors):
1) Lines 670672: The main assumption of the model is that nearsimultaneous binding of 2 calcium ions greatly increases the affinity of PIP2 binding. However, simultaneous binding of 2 Ca^{2+} ions must imply an even stronger cooperativity between the two Ca^{2+} binding events than the cooperativity between Ca^{2+} and PIP2 binding. While equilibrium properties are not affected by such an assumption, the timedependence of the model should depend on the precise sequence of binding events. Would the release latency predictions change if both Ca^{2+} binding events were resolved? Of course, the stochastic simulations for a more complicated model would be costly, and extra model parameters are not desired, but one trial analytic solution for a constant Ca^{2+} level is probably not hard. I am not insisting that any new simulations be included in the paper, but it would help to discuss briefly whether the assumption of two simultaneous Ca^{2+} binding affects the release latency estimation.
2) On page 10, the Authors describe their checks to ensure that the model parameter optimization procedure is not stuck at a local minimum. However, a more straightforward common check would involve repeating the fminsearch parameter optimization algorithm multiple times, while starting at different initial values of parameters. Have the Authors performed such a check? This would essentially quantify uncertainties in inferred parameter values in Table 2.
3) Apart from the nice connection between the proposed model and the model of Lou, Scheuss and Schneggenburger (2005), I think the presented model can also be viewed as a more detailed and biophysicallybased extension of the "excesscalcium binding site model" of S.D. Meriney and coworkers, which I would recommend citing. To my knowledge, there are 3 papers that make use of the latter model, in one form or another (but I would recommend that the Authors doublecheck these papers to see if all of them are relevant):
Dittrich M, Pattillo JM, King JD, Cho S, Stiles JR, Meriney SD (2013) An excesscalcium binding site model predicts neurotransmitter release at the neuromuscular junction. Biophys J 104: 275163
Ma J, Kelly L, Ingram J, Price TJ, Meriney SD, Dittrich M (2015) New insights into shortterm synaptic facilitation at the frog neuromuscular junction. J Neurophysiol 113: 7187
Luo F, Dittrich M, Cho S, Stiles JR, Meriney SD (2015) Transmitter release is evoked with low probability predominately by calcium flux through single channel openings at the frog neuromuscular junction. J Neurophysiol 113: 24809
Reviewer #2 (Recommendations for the authors):
1. The style and length of this paper resemble a thesis. The authors are encouraged to edit the manuscript to make it accessible to a wide eLife audience. Many figures have complicated subpanels that are difficult to understand, and the text is often so technical that it does not convey the essential point of each argument. This is a very nice model. The authors could do a better job of describing the important points, which they highlight most clearly at the beginning of the discussion.
2. The model depicted in Figure 1A suggests that Cabound C2B attaches to the vesicle, rather than the PM. However, Chang, Trimbuch and Rosenmund 2019 showed that the PIP2 binding attaches vesicles to the PM before an action potential. Do the authors want to claim that Cabound C2B attaches to vesicles? This is actually very important to the main point of the paper, because "crosslinks" are likely to occur in the absence of Ca, and Caindependent PIP2 binding might position the Cabinding pocket closer to negatively charged phospholipids, increasing the affinity of C2B for Ca.
3. In the discussion (Line 542) the authors note that their model is fundamentally different from the proposed preassembly of Syt rings at the base of a vesicle, and that this is big testable prediction. How could this be tested? And do Syts actually preassemble with slots in their model? Figure 3B shows the number of crosslinks that precede fusion, and it would be nice to see how many Syts are bound to PIP2 at rest for any given number of "slots".
4. Line 149: What (little) is known about the concentration of PIP2 at active zones? How does the concentration of PIP2 in rich patches compare to the idea that only 3 PIP2 molecules could be available beneath the space of one vesicle? I found the discussion of PIP2 concentration in Methods (Line 999) difficult to understand.
5. The model only accounts for Ca binding by the C2B domain, which as the authors note is the more important of the C2 domains. C2A is barely discussed in this paper. However, there are differences in fusion rates when the C2A domain is mutated to block Ca binding. How would the addition of C2A affect the model?
6. I like Figure 2 Sup2, where the authors probe the number of crosslinks prior to fusion. Interestingly, the speed of Ca ramping affects the number of crosslinks, and a significant number of events occur after 4 crosslinks only when Ca increases at the fastest rate (0.001 s). This is still much slower than the expected rate of Ca nanodomain formation (<100 us). Why did the authors explore only extremely slow Ca ramp rates? It seems that the ramping rates of 10 us to 10 ms would more accurately affect the condition of action potential induced Ca influx.
7. It is not immediately clear to me why the model with 6 slots never predicts fusion with 56 crosslinks. Why does the model fail to fit well if each crosslink contributes less to lowering the energy barrier, and thus 56 crosslinks are required to drive fusion?
8. If a single crosslink brings the vesicle closer to the PM, this will introduce a new form of allostery by increasing the effective concentration of PIP2 sensed by Syts. Is this reflected in the model?
9. In the model, PIP2 binding is required for Syt1driven fusion. As a result, only Syts with a "slot" can participate in fusion. However, while K325, 327A mutations that disrupt Caindependent PIP2 binding lead to less synchronous release, synchronous release can be restored using a pairedpulse stimulation paradigm where an initial action potential drives Cadependent membrane attachment of vesicles to the plasma membrane. This suggests that PIP2 binding is not required for "cross linking". How might the model change if Syts without slots could also contribute to fusion in your model?
10. The evaluation of mutants in Figure 5 should be better tied to actual biology. What mutation does "Ca^{2+} binding", and "Aon" mimic? No citations are provided for the many studies where Syt mutants were expressed in Syt1 KO neurons. This seems like a perfect place to test the role of PIP2 binding with more complicated Ca stimuli. Could the model be adapted to explore the effect of K/A mutations described by Chang, Trimbuch and Rosenmund?
Reviewer #3 (Recommendations for the authors):
Such kind of complex models tend to have various underlying assumptions which could significantly influence the conclusions. Three examples are provided:
1. The interaction of synaptotagmins with SNARES or other proteins could change the Ca^{2+} affinity but only PIP2 is considered as a potential interaction for changing the Ca^{2+} affinity. In other word, the premise of the paper that the Ca^{2+} affinity of synaptotagmin in vitro and in vivo (within the protein complex of the fusion machinery) is identical could be wrong. Therefore, the main conclusion of the strong cooperativity of Ca^{2+} and PIP2 could be wrong.
2. The assumption that each crosslinking synaptotagmin lowers the energy barrier for fusion by the same amount (E_syt) could be wrong. Assuming a positive or negative cooperativity in E_syt per synaptotagmin might impact the conclusions regarding the number of required synaptotagmins.
3. The distance between the individual synaptotagmin molecules and the nearest Ca^{2+} channels could differ significantly on the nmscale. The assumption of an identical local Ca^{2+} signal for all synaptotagmins could be wrong which could complicate the model predictions of the effect of reduced synaptotagmin copy numbers (Figure 4 and 5).
In conclusion, the study provides interesting and plausible possibilities of how synaptotagmins could mediate vesicle fusion but experimental validations are required to increase the amount of reliable and novel insights.
https://doi.org/10.7554/eLife.74810.sa1Author response
Essential revisions:
1) Lines 670672: The main assumption of the model is that nearsimultaneous binding of 2 calcium ions greatly increases the affinity of PIP2 binding. However, simultaneous binding of 2 Ca^{2+} ions must imply an even stronger cooperativity between the two Ca^{2+} binding events than the cooperativity between Ca^{2+} and PIP2 binding. While equilibrium properties are not affected by such an assumption, the timedependence of the model should depend on the precise sequence of binding events. Would the release latency predictions change if both Ca^{2+} binding events were resolved? Of course, the stochastic simulations for a more complicated model would be costly, and extra model parameters are not desired, but one trial analytic solution for a constant Ca^{2+} level is probably not hard. I am not insisting that any new simulations be included in the paper, but it would help to discuss briefly whether the assumption of two simultaneous Ca^{2+} binding affects the release latency estimation.
In our model, the binding of the two Ca^{2+} ions is described in one reaction, and the biological implications are connected to the state in which both Ca^{2+} ions have bound. This generates a simple interpretation regarding both the allosteric effect and the effect the C2B domains exert on the energy barrier for fusion:
1. The allosteric coupling between Ca^{2+} and PI(4,5)P_{2} is active when the Ca^{2+} binding sites are saturated (otherwise it is inactive).
2. Whenever a C2B domain binds both Ca^{2+} and PI(4,5)P_{2} the energy barrier for fusion is lowered.
If we split up the binding of Ca^{2+} into two consecutive reactions, we are faced with the problem of how to assign these two effects. For instance, is the allosteric effect “divided” between the Ca^{2+} bindings, or is it in effect only when both ions are bound? Does association of one Ca^{2+} together with PI(4,5)P_{2} have a partial effect on the energy barrier for fusion? How are the effects split up? As we do not know the answers to these questions, we circumvented these by assuming a simpler case in which both Ca^{2+} ions needed to be bound to exert these effects.
We agree with the reviewer that the exact mechanism is likely more complex. The coordination of the Ca^{2+} ions likely induces some conformational changes (e.g. membrane insertion) that stabilize the Ca^{2+} association and when this happens upon binding a single Ca^{2+} ion, the affinity for binding a second ion would likely be much higher. By simplifying the two binding reactions to one, we indeed indirectly assumed a high cooperativity between the binding events. We now added a statement of this to the methods section of the manuscript.
2) On page 10, the Authors describe their checks to ensure that the model parameter optimization procedure is not stuck at a local minimum. However, a more straightforward common check would involve repeating the fminsearch parameter optimization algorithm multiple times, while starting at different initial values of parameters. Have the Authors performed such a check? This would essentially quantify uncertainties in inferred parameter values in Table 2.
Repeating fminsearch with different initial parameters is a common way to confirm detection of the global minimum. Indeed, we had performed such an approach initially and have now also added a more systematic investigation. However, in our manuscript, we aimed to give more insight in the optimization space and therefore presented the analyses shown in figure 2F and 2G. These analyses show that our optimization problem has (within the parameter space that we explored) a global minimum at a cost of 139.45 and a local minimum with a cost of 107.
To address the reviewer’s comment we now fitted the model with 100 random initial values in which α, γ and F varied between 1 and 1000, [PI(4,5)P_{2}] between 1 and 100. The initial value for the added delay was set to 0.4 ms, as we approximately know its value. 40 out of 100 fits ended in the reported global minimum (cost = 139.45). When running the optimization procedure a second time with the fitted parameters as new initial values, a total of 60/100 fits ended in the reported minimum (139.45). During this fitting procedure, we also found a total of 18 fits with cost = 107.15 and a fitted fvalue of 989. Given that we observe a local minimum in figure 2H at f~1000, it is not surprising that multiple fits also converged to this minimum. The remaining 22 fits are likely not to have converged yet. Thus, running the fminsearch multiple times shows that the analysis performed in figure 2G and 2H, gives a good insight in the optimization problem.
3) In the discussion (Line 542) the authors note that their model is fundamentally different from the proposed preassembly of Syt rings at the base of a vesicle, and that this is big testable prediction. How could this be tested? And do Syts actually preassemble with slots in their model? Figure 3B shows the number of crosslinks that precede fusion, and it would be nice to see how many Syts are bound to PIP2 at rest for any given number of "slots".
The model of the preassembled syt rings (see for example Rothman et al., 2017) implies that rather specific quantities of synaptic proteins are needed to induce fusion. Such a model would be highly sensitive to reducing the number of participating syts. Transmitter release would immediately break down if the number of syts were reduced such that no rings could preassemble. This contrasts our model where titration of syt copy number per SV yields a gradual effect on release and where a high number of syt copies is not needed for SV fusion per se but ensures fast and reliable neurotransmission by increasing the speed of the response to a stimulus (Figure 5). This difference suggests a way to test these hypotheses experimentally. We extended our Discussion section on this.
Whether syts preassemble with slots/PI(4,5)P_{2} in our model is a very interesting question. We therefore added the results on the steady state distribution as an additional figure (figure 3, Figure 3 —figure supplement 1) to our manuscript. Besides looking into the steady state distribution, we also investigated whether those vesicles that have multiple syts associated to PI(4,5)P_{2} have a higher probability of contributing to fusion. Below, we describe these results.
The histogram in Author response image 1 illustrates the steady state for each choice of M_{slots} (see also Figure 3A, and figure 3 —figure supplement 1A for steady state distributions of M_{slots} = 3 and M_{slots} = 6 respectively). For each choice of M_{slots}, we used the corresponding best fit parameters, which makes parameters vary between choices of M_{slots.} The bars show the percentage of SVs having bound any number of PI(4,5)P_{2} molecules at steady state without forming any crosslinks. Formation of a crosslink at steady state happens with probability between 0.02 % and 0.0001 % depending on the number of slots.
This shows that for M_{slots}≥3 a considerable number of SVs preassembles with PI(4,5)P_{2} prior to the stimulus. In the following plots we investigated the consequence of steady state PI(4,5)P_{2} binding for the best fit situation, that is for M_{slots}=3.An SV with one or more syts having bound PI(4,5)P_{2} at steady state naturally can have an “advantage” in release probability. This is illustrated in Author response image 2, where the release probability of SVs with a certain amount of PI(4,5)P_{2} bound at steady state is plotted against time after a Ca^{2+} flash of 50 µM see also Figure 3C (M_{slots} = 3), and Figure 3—figure supplement 1C (M_{slots}=6).
Combining the steady state amounts with the release probabilities, we can compute the cumulative number of SV fusions over time with an RRP size of 4000 SVs (Figure 3B).Relatively many of the fastest SV fusions are a result of syts preassembling with PI(4,5)P_{2}. Restricting the plot to the first 0.05 ms, we see that preassembling with PI(4,5)P_{2} really affects the initial response to this stimulus, which is highly relevant for the obtained release latency (defined as the time between stimulus and the 5^{th} SV fusion event).
Since the affinities for Ca^{2+} and PI(4,5)P_{2} are fixed to experimental values, the PI(4,5)P_{2} concentration is the only fitted parameter that can actually change the PI(4,5)P_{2} binding status at steady state. This provides an explanation for the large impact of the PI(4,5)P_{2} concentration on the release latencies (as seen in Figure 2 —figure supplement 5) and justifies the fitting of both the PI(4,5)P_{2} concentration and the PI(4,5)P_{2} binding rate, as these parameters affect the model in very different ways.
The above analyses are performed on the level of the vesicles. We also wanted to investigate the relevance of having bound to PI(4,5)P_{2} at steady state at the level of individual syts.
For this, like in some parts of the manuscript, we used the Gillespie algorithm. We ran the algorithm with 100 vesicles 200 repetitions for the Calcium steps to 1, 5, 10, 50, and 100 µM. As we were interested in investigating the probability of being involved in fusion (engaging in Ca^{2+}/PI(4,5)P_{2} dualbinding) for those syts that have prebound to PI(4,5)P_{2} at steady state, we only ran the algorithm until all syts that were bound to PI(4,5)P_{2} at steady state had either fused or unbound PI(4,5)P_{2}. In Author response table 1 you can find the proportion of syts prebound to PI(4,5)P_{2} at steady state that is involved in fusion for the different Calcium concentrations analysed. In line with the above results, at a calcium concentration of 50 µM, the vast majority of the syts that have bound PI(4,5)P_{2} at steady state are also actively contributing to fusion by dual binding Ca^{2+}/PI(4,5)P_{2}. The proportion of prebounds syts involved in fusion is highly dependent on the Calcium concentration. At low Calcium concentrations (1 and 5 µM) having PI(4,5)P_{2} bound at steady state is not a predictor for being involved infusion.
4) I like Figure 2 Sup2, where the authors probe the number of crosslinks prior to fusion. Interestingly, the speed of Ca ramping affects the number of crosslinks, and a significant number of events occur after 4 crosslinks only when Ca increases at the fastest rate (0.001 s). This is still much slower than the expected rate of Ca nanodomain formation (<100 us). Why did the authors explore only extremely slow Ca ramp rates? It seems that the ramping rates of 10 us to 10 ms would more accurately affect the condition of action potential induced Ca influx.
This relies on a misunderstanding. We did explore very fast (i.e. instantaneous behaviour) in the left hand column and then explored slower ramps in the right hand column. At the fastest rise time (0.001 s), convergence was already reached, as the average of approx. 3.4 crosslinks before fusion was also seen with 100 µM step Ca^{2+}. We have now tried to make this point clearer by adding a few faster rise time points, making the xaxis logarithmic as well as drawing a horizontal line that visualises the convergence.
5) It is not immediately clear to me why the model with 6 slots never predicts fusion with 56 crosslinks. Why does the model fail to fit well if each crosslink contributes less to lowering the energy barrier, and thus 56 crosslinks are required to drive fusion?
Indeed, it is striking that we consistently see 34 dual bindings formed before fusions – even when allowing up to 6 dualbounds. Fitting the model with different choices of M_{slots}≥3 yielded consistent estimation of the contribution to the energy barrier for fusion by each syt engaging in Ca^{2+}/PI(4,5)P_{2} dual binding (Figure 2F), even when probing many different initial parameter values in the fitting routine. As the reviewer notes, if each individual dualbinding syt contributed less to the energy barrier for fusion, more C2B domains would be required to engage simultaneously to achieve the fast rates of SV fusion. However, engaging more domains also would require the binding of more Ca^{2+} ions which could alter the Ca^{2+} dependence of the release rates. To investigate this, we ran the fitting procedure forcing the model with 6 slots to induce fusion with fastest rate only after 6 dualbounds have been formed. This we did by fixing the contribution of each individual dualbinding syt to the energy barrier (factor f) during the fitting procedure and only fitted the four remaining parameters for a model with M_{slots} = 6. We computed the fixed f value based on the f value obtained by fitting all five model parameters, which equalled 163.5 for the six slots condition (M_{slots} = 6). The f value that, under the same conditions, would lead to most fusion events to occur once 5/6 dualbounds have formed can be computed by: ${\left(163.5\right)}^{\frac{3}{6}}=12.79$ (see Figure 2 —figure supplement 4 for confirmation). The best fit that could be obtained using f = 12.79 had a cost value of 8.78, which is much larger than the best fit obtained by fitting with f as a free parameter (cost = 126.5). Figure 2 —figure supplement 4 (panel A) show the best fit release latencies and peak release rates when forcing the six syts engaging in Ca^{2+}/PI(4,5)P_{2} dual binding during fast fusion. The 95% prediction interval captures less of the data points for both the release latency and peak release rates condition compared to the original fit. Furthermore, as expected, the model predicted a calciumdependency of the peak release rates which was too steep compared to the experimental data. Panel B shows that with these new settings, most fusion events occur when 6 crosslinks have formed. We added this figure as an additional supplement to our manuscript (Figure 2 —figure supplement 4)
The requirement for three crosslinking syts for fusion likely relates to the assumption that each C2B domain binds two Ca^{2+} ions. We also tested whether the condition with six slots and the binding of a single Ca^{2+} ion to each synaptotagmin could result in model predictions corresponding to the experimental data (Ca^{2+} affinity and allostericity value adjusted such that they correspond to binding of a single Ca^{2+} ion), and show fusion from a condition were six synaptotagmins are binding both Ca^{2+} and PI(4,5)P_{2}. Also with these settings (using K_{D,Ca}^{2+} = 221 µM, and A = 0.0015), the model was not fitting to the experimental data (see Author response image 3 and Author response table 2).
6) If a single crosslink brings the vesicle closer to the PM, this will introduce a new form of allostery by increasing the effective concentration of PIP2 sensed by Syts. Is this reflected in the model?
This is a very good point and we have added a new Figure (New Figure 7) to indicate the strong effect changes of the effective PI(4,5)P_{2} concentration has on synaptic responses (see point below). Currently our model does not include any heterogeneity in the effective PI(4,5)P_{2} concentration individual syts experience. This could be the case if all RRP vesicles are similarly docked. Indeed, our model predicts that most RRP SVs preassociate PI(4,5)P_{2} (New Figure 3). However, it remains possible that binding of additional SV syts to PI(4,5)P_{2} further affects the effective PI(4,5)P_{2} concentration sensed by the remaining syts. We currently do not see how we might quantitatively estimate this effect. As a matter of fact, the effect may be bidirectional for the remaining syts on SV (e.g. while the effective PI(4,5)P_{2} concentration may increase for syts facing same PI(4,5)P_{2}patch on the plasma membrane, it may decrease for ones facing away from the PI(4,5)P_{2} patch). We currently do not know how the size of these effects may be estimated (this would require taking into account geometrical information, diffusion of the syt transmembrane anchor, protein lengths and flexibility, and the organization of PI(4,5)P_{2} in the membrane). We therefore kept our model as simple as possible in this regard, by assuming that all syst see the same PI(4,5)P_{2} concentration. Thus, the best fit PI(4,5)P_{2} concentration may reflect the average conditions for each syt. We now explain this clearer in the Results and touch on the above limitation in the Discussion.
7) In the model, PIP2 binding is required for Syt1driven fusion. As a result, only Syts with a "slot" can participate in fusion. However, while K325, 327A mutations that disrupt Caindependent PIP2 binding lead to less synchronous release, synchronous release can be restored using a pairedpulse stimulation paradigm where an initial action potential drives Cadependent membrane attachment of vesicles to the plasma membrane. This suggests that PIP2 binding is not required for "cross linking". How might the model change if Syts without slots could also contribute to fusion in your model?
To address this point, we performed additional simulations and have now added a new figure (Figure 7) to the manuscript where we explore the consequences of activity dependent SV repositioning. This point also nicely links to the one above regarding the effective PI(4,5)P_{2} concentration. As pointed out above, our model features one PI(4,5)P_{2} concentration that determines the likelihood that syts associate. This PI(4,5)P_{2} concentration encompasses both the density of PI(4,5)P_{2} in the membrane and the accessibility of syt to PI(4,5)P_{2}. The latter will depend on how close the vesicle is docked to the plasma membrane. In their paper, Chang, Timbruch and Rosenmund (Chang et al., 2018) show that mutation of basic residues of the C2B domain (lysines also implicated in PI(4,5)P_{2} binding or the aspartates implicated in membrane binding) leads to loss of SV docking. Docking is restored upon APstimulation of mutant synapses and probing synaptic responses with a second AP 10 ms after the first leads to a large facilitation of synaptic responses. We here simulated such a condition by implementing a timedependent increase in the effective PI(4,5)P_{2} concentration in such mutants. The assumption here is that when SVs are distant from the PM, the effective PI(4,5)P_{2} level available for syts are low (“the membrane is more difficult to reach”) but this situation normalizes upon SV translocation to the PM (Figure 7A). Accordingly, responses to the first AP are strongly reduced, but normalize for a second AP, causing strong shortterm facilitation (Figure 7BD). While our model clearly is a strong simplification (as stated above, the effective PI(4,5)P_{2} concentration in reality is a property relevant for each individual syt and this situation becomes more relevant here due to the more inhomogeneous distribution of RRP SVs) our model can provide a simple interpretation of the phenotypic cause of this type of mutations. We thank the reviewer for giving us the opportunity to explore this interesting point.
8) The evaluation of mutants in Figure 5 should be better tied to actual biology. What mutation does "Ca^{2+} binding", and "Aon" mimic? No citations are provided for the many studies where Syt mutants were expressed in Syt1 KO neurons. This seems like a perfect place to test the role of PIP2 binding with more complicated Ca stimuli. Could the model be adapted to explore the effect of K/A mutations described by Chang, Trimbuch and Rosenmund?
We agree that this section of the manuscript was not presented well. We have now rewritten the text dealing with the putative mutations that interfere with Ca^{2+} binding and the allosteric coupling between the Ca^{2+} and PI(4,5)P_{2} binding sites. We now specifically name and reference the mutations that our hypothetical model parameter changes may relate to. We have furthermore reduced the complexity of this figure (removed superfluous panels/supplementary items).
Following the reviewers’ suggestions (see point above) we have now also added an exploration of putative effects caused by mutations that interfere with plasma membrane association. In this context we also used the opportunity to provide further analysis on repetitive AP stimuli. The new results are added in a new section of the main text and one new main figure (Figure 7). Our simulations show that activitydependent changes in the effective PI(4,5)P_{2} concentration (caused by SV repositioning) will have prominent effects on shortterm synaptic plasticity. We thank the reviewer for this suggestion and for giving us the opportunity to provide these further insights.
9) Recent experimental and modeling work with Synaptotagmin and SNAREs has been published (Wu et al., 2021; eLife; https://elifesciences.org/articles/68215). These authors say: "To test whether Syt1 affected fusion pores in this system, we coreconstituted ~4 copies of recombinant fulllength Syt1 together with ~4 copies of VAMP2 (per disc face) into large nanodiscs called nanolipoprotein particles…". The work seems particularly relevant to your paper. Also please take a look at another paper from this group that suggests fast exocytosis requires up to 15 SNAREs and Syt complexes (Wu et al., eLife, 2017; https://elifesciences.org/articles/22964). Please read and discuss in your paper these recent modeling studies, which seem to suggest that large numbers of Syt's and SNAREs are needed for fast exocytosis at synapses.
We would like to thank the reviewer for pointing out these relevant papers. We have added the points to the discussion. We especially discuss that synaptotagmin was found to only change the fusion kinetics when PI(4,5)P_{2} was present and we included this finding in our manuscript to give a better explanation for why in our model only if both Ca^{2+} and PI(4,5)P_{2} are bound to the C2B domain the fusion barrier for fusion is reduced.
Reviewer #1 (Recommendations for the authors):
1) Lines 670672: The main assumption of the model is that nearsimultaneous binding of 2 calcium ions greatly increases the affinity of PIP2 binding. However, simultaneous binding of 2 Ca^{2+} ions must imply an even stronger cooperativity between the two Ca^{2+} binding events than the cooperativity between Ca^{2+} and PIP2 binding. While equilibrium properties are not affected by such an assumption, the timedependence of the model should depend on the precise sequence of binding events. Would the release latency predictions change if both Ca^{2+} binding events were resolved? Of course, the stochastic simulations for a more complicated model would be costly, and extra model parameters are not desired, but one trial analytic solution for a constant Ca^{2+} level is probably not hard. I am not insisting that any new simulations be included in the paper, but it would help to discuss briefly whether the assumption of two simultaneous Ca^{2+} binding affects the release latency estimation.
Answered in the list of essential revisions above.
2) On page 10, the Authors describe their checks to ensure that the model parameter optimization procedure is not stuck at a local minimum. However, a more straightforward common check would involve repeating the fminsearch parameter optimization algorithm multiple times, while starting at different initial values of parameters. Have the Authors performed such a check? This would essentially quantify uncertainties in inferred parameter values in Table 2.
Answered in the list of essential revisions above.
3) Apart from the nice connection between the proposed model and the model of Lou, Scheuss and Schneggenburger (2005), I think the presented model can also be viewed as a more detailed and biophysicallybased extension of the "excesscalcium binding site model" of S.D. Meriney and coworkers, which I would recommend citing. To my knowledge, there are 3 papers that make use of the latter model, in one form or another (but I would recommend that the Authors doublecheck these papers to see if all of them are relevant):
Dittrich M, Pattillo JM, King JD, Cho S, Stiles JR, Meriney SD (2013) An excesscalcium binding site model predicts neurotransmitter release at the neuromuscular junction. Biophys J 104: 275163
Ma J, Kelly L, Ingram J, Price TJ, Meriney SD, Dittrich M (2015) New insights into shortterm synaptic facilitation at the frog neuromuscular junction. J Neurophysiol 113: 7187
Luo F, Dittrich M, Cho S, Stiles JR, Meriney SD (2015) Transmitter release is evoked with low probability predominately by calcium flux through single channel openings at the frog neuromuscular junction. J Neurophysiol 113: 24809
We agree with the reviewer the excess calciumbinding site model has many similarities to our model. In our discussion, we now include a discussion of the original article. The followup papers of the model focus on facilitation and the origin of the Ca^{2+} ions that induce fusion, which we did not consider as relevant to our study.
Reviewer #2 (Recommendations for the authors):
1. The style and length of this paper resemble a thesis. The authors are encouraged to edit the manuscript to make it accessible to a wide eLife audience. Many figures have complicated subpanels that are difficult to understand, and the text is often so technical that it does not convey the essential point of each argument. This is a very nice model. The authors could do a better job of describing the important points, which they highlight most clearly at the beginning of the discussion.
We appreciate this suggestion and have extensively rewritten and simplified the text. We now put a stronger emphasis on the biological relevance and have moved more technical sections to the methods. We have also worked to simplify our figures.
2. The model depicted in Figure 1A suggests that Cabound C2B attaches to the vesicle, rather than the PM. However, Chang, Trimbuch and Rosenmund 2019 showed that the PIP2 binding attaches vesicles to the PM before an action potential. Do the authors want to claim that Cabound C2B attaches to vesicles? This is actually very important to the main point of the paper, because "crosslinks" are likely to occur in the absence of Ca, and Caindependent PIP2 binding might position the Cabinding pocket closer to negatively charged phospholipids, increasing the affinity of C2B for Ca.
We realized that our term “crosslink” and the illustration of our model did give the impression that the C2B domain must bind both the plasma membrane and the vesicular membrane to exert its action. However, this is only one possible mode of action and in fact our model does not rely on this assumption at all. We now tried to clarify this by directly pointing out that it is currently not known how the C2B exerts its effect on the fusion barrier and name several possible mechanisms (induction of curvature, switch in the local electrostatic environment, direct or indirect regulation of SNARE complex assembly and distance regulation). To clarify this have removed the term “crosslink” from our manuscript. We have also included and additional supplementary figure (Figure 1 —figure supplement 1) for another alternative that the C2B only associates to the plasma membrane.
3. In the discussion (Line 542) the authors note that their model is fundamentally different from the proposed preassembly of Syt rings at the base of a vesicle, and that this is big testable prediction. How could this be tested? And do Syts actually preassemble with slots in their model? Figure 3B shows the number of crosslinks that precede fusion, and it would be nice to see how many Syts are bound to PIP2 at rest for any given number of "slots".
Answered in the list of essential revisions above.
4. Line 149: What (little) is known about the concentration of PIP2 at active zones? How does the concentration of PIP2 in rich patches compare to the idea that only 3 PIP2 molecules could be available beneath the space of one vesicle? I found the discussion of PIP2 concentration in Methods (Line 999) difficult to understand.
Currently little is known about the exact density of PI(4,5)P_{2} at release sites. Indeed, we know that PI(4,5)P_{2} forms clusters on membranes, which contain a high density of PI(4,5)P_{2} (Honigmann et al., 2013; Milosevic et al., 2005; van den Bogaart et al., 2011a). Clustering of PI(4,5)P_{2} on the plasma membrane is one of the reasons why we implemented slots in our model, but there are also several other explanations why the number of syts simultaneously binding to PI(4,5)P_{2} is limited. The number of slots in our model, which is set to three, does not correspond to the number of PI(4,5)P_{2} molecules beneath the space of one vesicle. The slots in our model set a limit in the number of Syts that can bind to PI(4,5)P_{2} simultaneously. We assume that each slot contains a certain PI(4,5)P_{2} concentration, which allowed us to directly implement binding affinities determined in the in vitro experiment by van den Bogaart et al., to constrain our model. Importantly, as for all models with three or more slots, we observe that most fusion events occur when three dualbounds are formed, the number of slots in our model will not impact the main conclusions we draw. We now extend on the discussion of these concepts and the simplifications that were necessary to implement the model.
The section on the interpretation of the PI(4,5)P_{2} concentration in Methods was indeed quite technical. The main point we wanted to make in this section is that our model leads to an estimation of the concentration of PI(4,5)P_{2}, whereas an estimate of the density of PI(4,5)P_{2} in the membrane would be more relevant, as PI(4,5)P_{2} is not dissolved in the cytosol. However, to compute the density from the estimated concentration we would need to include many assumptions. We removed this section from Methods and implemented the points in the Discussion.
5. The model only accounts for Ca binding by the C2B domain, which as the authors note is the more important of the C2 domains. C2A is barely discussed in this paper. However, there are differences in fusion rates when the C2A domain is mutated to block Ca binding. How would the addition of C2A affect the model?
We agree with the reviewer that we should discuss the C2A domain of Synaptotagmin more in our paper. We therefore included an additional section in the discussion.
The number of Ca^{2+} ions needed to bind to Synaptotagmin to induce fusion is the only parameter in our model that is based specifically on the C2B domain. The used Ca^{2+} and PI(4,5)P_{2} affinities are obtained using a construct containing both the C2A and C2B domain. However as mutating the Ca^{2+} binding sites of the C2A domain did not significantly affect the measured affinities, those affinity values can be attributed to the C2B domain (van den Bogaart, et al., 2012). Moreover, as we fitted our model to experimental data which is obtained in the presence of both the C2A and C2B domain of Synaptotagmin, our fitted model parameters likely also include some properties of the C2A domain. Therefore, although we focus on the C2B domain of Synaptotagmin, properties of the C2A domain are most likely indirectly included in our model simulations. For instance, quantification of the Ca^{2+} sensitivity of dense core vesicle fusion revealed a change upon mutation of the C2A domain binding sites (Soerensen et al., 2003), indicating that Ca^{2+} binding to the C2A domain may affect the affinity of the C2B domain. A similar effect in our model is attributed to a change in the PI(4,5)P_{2} concentration. Thus, a contribution of the C2A domain may be indirectly represented by our determined parameters. We agree with the reviewer that the initial version of the manuscript contained too little discussion on the C2A domain. We therefore now wrote a section on the C2A domain in the discussion.
An explicit simulation of the role of the C2A domain would require the inclusion of three more Ca^{2+} binding pockets with other, unknown Ca^{2+} affinities. And it might require additional assumptions on how the energy barrier for synaptic vesicle fusion is affected if Ca^{2+} ions bind to the C2A domain. While this may indeed be relevant for synaptic function, our simpler model focussing on the C2B domain only was sufficient to account for the experimental data we evaluated in this study. We thus chose not to explicitly include Ca^{2+} binding to the C2A domain in the context of this model.
6. I like Figure 2 Sup2, where the authors probe the number of crosslinks prior to fusion. Interestingly, the speed of Ca ramping affects the number of crosslinks, and a significant number of events occur after 4 crosslinks only when Ca increases at the fastest rate (0.001 s). This is still much slower than the expected rate of Ca nanodomain formation (<100 us). Why did the authors explore only extremely slow Ca ramp rates? It seems that the ramping rates of 10 us to 10 ms would more accurately affect the condition of action potential induced Ca influx.
Answered in the list of essential revisions above.
7. It is not immediately clear to me why the model with 6 slots never predicts fusion with 56 crosslinks. Why does the model fail to fit well if each crosslink contributes less to lowering the energy barrier, and thus 56 crosslinks are required to drive fusion?
Answered in the list of essential revisions above.
8. If a single crosslink brings the vesicle closer to the PM, this will introduce a new form of allostery by increasing the effective concentration of PIP2 sensed by Syts. Is this reflected in the model?
Answered in the list of essential revisions above.
9. In the model, PIP2 binding is required for Syt1driven fusion. As a result, only Syts with a "slot" can participate in fusion. However, while K325, 327A mutations that disrupt Caindependent PIP2 binding lead to less synchronous release, synchronous release can be restored using a pairedpulse stimulation paradigm where an initial action potential drives Cadependent membrane attachment of vesicles to the plasma membrane. This suggests that PIP2 binding is not required for "cross linking". How might the model change if Syts without slots could also contribute to fusion in your model?
Answered in the list of essential revisions above.
10. The evaluation of mutants in Figure 5 should be better tied to actual biology. What mutation does "Ca^{2+} binding", and "Aon" mimic? No citations are provided for the many studies where Syt mutants were expressed in Syt1 KO neurons. This seems like a perfect place to test the role of PIP2 binding with more complicated Ca stimuli. Could the model be adapted to explore the effect of K/A mutations described by Chang, Trimbuch and Rosenmund?
Answered in the list of essential revisions above.
Reviewer #3 (Recommendations for the authors):
Such kind of complex models tend to have various underlying assumptions which could significantly influence the conclusions. Three examples are provided:
All models are based on assumptions and rested on observations. In all assumptions we made we aimed to match currently available experimental evidence as closely as possible. Additionally, we tried to keep the model as simple as possible to limit the number of assumptions needed to be made. To add strength to the model predictions performance of experiments would be of great value. However, some of the model predictions, such as the importance of the allosteric interaction for the function of Synaptotagmin, currently cannot be tested experimentally. This is often the motivation of deriving a model which can make testable predictions. We therefore believe that our modelling approach provides unique insights in the working mechanism of synaptotagmins.
1. The interaction of synaptotagmins with SNARES or other proteins could change the Ca^{2+} affinity but only PIP2 is considered as a potential interaction for changing the Ca^{2+} affinity. In other word, the premise of the paper that the Ca^{2+} affinity of synaptotagmin in vitro and in vivo (within the protein complex of the fusion machinery) is identical could be wrong. Therefore, the main conclusion of the strong cooperativity of Ca^{2+} and PIP2 could be wrong.
We respectfully point out that models are not typically “judged” on whether all assumptions and predictions are going turn out to be “right” in the future, rather by their use to make a synthesis of current knowledge and to make predictions inaccessible to measurements. Clearly, with future insights also our model will need to be extended.
We agree with the reviewer that interactions with the SNARE complex could also potentially increase the Ca^{2+} affinity of synaptotagmin by binding to its polybasic sequence via a similar mechanisms as PI(4,5)P_{2} binding. However, such a situation would still be consistent with our model and the strong allosteric mechanism between the Ca^{2+} and PI(4,5)P_{2} and/or SNARE binding sites essential for the protein’s function (as we point out). Therefore, this point is not at odds with our model and the Ca^{2+} affinity in vitro and vivo would match, simply the interpretation of the molecular interaction responsible for the Ca^{2+} affinity increase differs.
2. The assumption that each crosslinking synaptotagmin lowers the energy barrier for fusion by the same amount (E_syt) could be wrong. Assuming a positive or negative cooperativity in E_syt per synaptotagmin might impact the conclusions regarding the number of required synaptotagmins.
It is a possibility that the additive contribution of individual syts to the energy barrier is nonlinear. However, even if this were the case, there is currently no experimental measurement that would allow to constrain such a complex scenario, even the direction of the cooperativity (negative or positive) cannot be predicted. Implementation of such further mechanisms to the model would go against this reviewer’s primary concern as it would require additional assumptions and unknown parameters. We avoid this uncertainty by investigating the simplest scenario (independent and additive effects).
In case of a positive cooperativity, i.e. each dualbound syt having an increasing effect on the effect on the height of the energy barrier for synaptic vesicle fusion, our model would still predict the involvement of at least two synaptotagmins in fusion, as we could not fit our model with only 1 slot (Figure 2). With a negative cooperativity, i.e. each dualbound synaptotagmin has a reducing effect on the reduction of the height of the energy barrier, the number of synaptotagmins actively involved in fusion could be larger than three. Similar to the points raised above, such specialisations would infer slight modifications of the model, but would not lead to a different model per se.
3. The distance between the individual synaptotagmin molecules and the nearest Ca^{2+} channels could differ significantly on the nmscale. The assumption of an identical local Ca^{2+} signal for all synaptotagmins could be wrong which could complicate the model predictions of the effect of reduced synaptotagmin copy numbers (Figure 4 and 5).
The data the model was fit to were obtained in Ca^{2+} uncaging experiments. In these experiments, the distance to voltage gated Ca^{2+} channels does not matter because the stimulus is spatially homogenous. Also, to simulate the action potential evoked EPSCs, we used a previously published calcium transient that was calculated from calciumuncaging data and measured EPSC waveforms, while assuming that all vesicles experience the same Ca^{2+}signal (Wang, Neher, Taschenberger 2008). This can be viewed as the ‘average’ Ca^{2+} signal RRP vesicles experience upon arrival of an action potential. We agree that there can be interesting effects of combining a heterogeneous vesiclesynaptotagmin number with heterogeneous vesiclechannel distances, but this will again require additional unknown parameters.
In conclusion, the study provides interesting and plausible possibilities of how synaptotagmins could mediate vesicle fusion but experimental validations are required to increase the amount of reliable and novel insights.
We thank the reviewer for the encouraging remarks and in our manuscript now point out several experiments that could be used to test different aspects of the model.
https://doi.org/10.7554/eLife.74810.sa2Article and author information
Author details
Funding
Novo Nordisk Fonden (NNF19OC0056047)
 Alexander M Walter
Novo Nordisk Fonden (NNF20OC0062958)
 Susanne Ditlevsen
Novo Nordisk Fonden (NNF19OC0058298)
 Jakob B Sørensen
Lundbeckfonden (R2772018802)
 Jakob B Sørensen
Independent Research Fund Denmark (802000228A)
 Jakob B Sørensen
Deutsche Forschungsgemeinschaft (278001972)
 Alexander M Walter
Deutsche Forschungsgemeinschaft (261020751)
 Alexander M Walter
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Ralf Schneggenburger and Holger Taschenberger for providing the experimental Ca^{2+} uncaging dataset and APinduced Ca^{2+} transient.
Senior Editor
 Richard W Aldrich, The University of Texas at Austin, United States
Reviewing Editor
 Henrique von Gersdorff, Oregon Health and Science University, United States
Reviewer
 Victor Matveev, New Jersey Institute of Technology, United States
Publication history
 Received: October 18, 2021
 Preprint posted: October 23, 2021 (view preprint)
 Accepted: August 4, 2022
 Accepted Manuscript published: August 5, 2022 (version 1)
 Version of Record published: September 20, 2022 (version 2)
Copyright
© 2022, Kobbersmed, Berns et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 885
 Page views

 354
 Downloads

 1
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Computational and Systems Biology
Computational models starting from large ensembles of evolutionarily related protein sequences capture a representation of protein families and learn constraints associated to protein structure and function. They thus open the possibility for generating novel sequences belonging to protein families. Protein language models trained on multiple sequence alignments, such as MSA Transformer, are highly attractive candidates to this end. We propose and test an iterative method that directly employs the masked language modeling objective to generate sequences using MSA Transformer. We demonstrate that the resulting sequences score as well as natural sequences, for homology, coevolution, and structurebased measures. For large protein families, our synthetic sequences have similar or better properties compared to sequences generated by Potts models, including experimentally validated ones. Moreover, for small protein families, our generation method based on MSA Transformer outperforms Potts models. Our method also more accurately reproduces the higherorder statistics and the distribution of sequences in sequence space of natural data than Potts models. MSA Transformer is thus a strong candidate for protein sequence generation and protein design.

 Cancer Biology
 Computational and Systems Biology
Lung squamous cell carcinoma (LUSC) is a type of lung cancer with a dismal prognosis that lacks adequate therapies and actionable targets. This disease is characterized by a sequence of low and highgrade preinvasive stages with increasing probability of malignant progression. Increasing our knowledge about the biology of these premalignant lesions (PMLs) is necessary to design new methods of early detection and prevention, and to identify the molecular processes that are key for malignant progression. To facilitate this research, we have designed XTABLE (Exploring Transcriptomes of Bronchial Lesions), an opensource application that integrates the most extensive transcriptomic databases of PMLs published so far. With this tool, users can stratify samples using multiple parameters and interrogate PML biology in multiple manners, such as two and multiplegroup comparisons, interrogation of genes of interests, and transcriptional signatures. Using XTABLE, we have carried out a comparative study of the potential role of chromosomal instability scores as biomarkers of PML progression and mapped the onset of the most relevant LUSC pathways to the sequence of LUSC developmental stages. XTABLE will critically facilitate new research for the identification of early detection biomarkers and acquire a better understanding of the LUSC precancerous stages.