Abstract
Chemical synaptic transmission relies on the Ca^{2+}induced fusion of transmitterladen vesicles whose coupling distance to Ca^{2+} channels determines synaptic release probability and shortterm plasticity, the facilitation or depression of repetitive responses. Here, using electron and superresolution microscopy at the Drosophila neuromuscular junction we quantitatively map vesicle:Ca^{2+} channel coupling distances. These are very heterogeneous, resulting in a broad spectrum of vesicular release probabilities within synapses. Stochastic simulations of transmitter release from vesicles placed according to this distribution revealed strong constraints on shortterm plasticity; particularly facilitation was difficult to achieve. We show that postulated facilitation mechanisms operating via activitydependent changes of vesicular release probability (e.g. by a facilitation fusion sensor) generate too little facilitation and too much variance. In contrast, Ca^{2+}dependent mechanisms rapidly increasing the number of releasable vesicles reliably reproduce shortterm plasticity and variance of synaptic responses. We propose activitydependent inhibition of vesicle unpriming or release site activation as novel facilitation mechanisms.
eLife digest
Cells in the nervous system of all animals communicate by releasing and sensing chemicals at contact points named synapses. The ‘talking’ (or presynaptic) cell stores the chemicals close to the synapse, in small spheres called vesicles. When the cell is activated, calcium ions flow in and interact with the releaseready vesicles, which then spill the chemicals into the synapse. In turn, the ‘listening’ (or postsynaptic) cell can detect the chemicals and react accordingly.
When the presynaptic cell is activated many times in a short period, it can release a greater quantity of chemicals, allowing a bigger reaction in the postsynaptic cell. This phenomenon is known as facilitation, but it is still unclear how exactly it can take place. This is especially the case when many of the vesicles are not ready to respond, for example when they are too far from where calcium flows into the cell. Computer simulations have been created to model facilitation but they have assumed that all vesicles are placed at the same distance to the calcium entry point: Kobbersmed et al. now provide evidence that this assumption is incorrect.
Two highresolution imaging techniques were used to measure the actual distances between the vesicles and the calcium source in the presynaptic cells of fruit flies: this showed that these distances are quite variable – some vesicles sit much closer to the source than others.
This information was then used to create a new computer model to simulate facilitation. The results from this computing work led Kobbersmed et al. to suggest that facilitation may take place because a calciumbased mechanism in the cell increases the number of vesicles ready to release their chemicals.
This new model may help researchers to better understand how the cells in the nervous system work. Ultimately, this can guide experiments to investigate what happens when information processing at synapses breaks down, for example in diseases such as epilepsy.
Introduction
At chemical synapses, neurotransmitters (NTs) are released from presynaptic neurons and subsequently activate postsynaptic receptors to transfer information. At the presynapse, incoming action potentials (APs) trigger the opening of voltage gated Ca^{2+} channels, leading to Ca^{2+} influx. This local Ca^{2+} signal induces the rapid fusion of NTcontaining synaptic vesicles (SVs) at active zones (AZs) (Südhof, 2012). In preparation for fusion, SVs localize (dock) to the AZ plasma membrane and undergo functional maturation (priming) into a readily releasable pool (RRP) (Kaeser and Regehr, 2017; Verhage and Sørensen, 2008). These reactions are mediated by an evolutionarily highly conserved machinery. The SV protein VAMP2/Synaptobrevin and the plasma membrane proteins Syntaxin1 and SNAP25 are essential for docking and priming and the assembly of these proteins into the ternary SNARE complex provides the energy for SV fusion (Jahn and Fasshauer, 2012). The SNARE interacting proteins (M)Unc18s and (M)Unc13s (where ‘M’ indicates mammalian) are also essential for SV docking, priming and NT release (Rizo and Südhof, 2012; Südhof and Rothman, 2009), while Ca^{2+} triggering of SV fusion depends on vesicular Ca^{2+} sensors of the Synaptotagmin family (Littleton and Bellen, 1995; Südhof, 2013; Walter et al., 2011; Yoshihara et al., 2003). Cooperative binding of multiple Ca^{2+} ions to the SV fusion machinery increases the probability of SV fusion (pV_{r}) in a nonlinear manner (Bollmann et al., 2000; Dodge and Rahamimoff, 1967; Schneggenburger and Neher, 2000).
A distinguishing feature of synapses is their activity profile upon repeated AP activation, where responses deviate between successive stimuli, resulting in either shortterm facilitation (STF) or shortterm depression (STD). This shortterm plasticity (STP) fulfils essential temporal computational tasks (Abbott and Regehr, 2004). Postsynaptic STP mechanisms can involve altered responsiveness of receptors to NT binding, while presynaptic mechanisms can involve alterations in Ca^{2+} signalling and –sensitivity of SV fusion (von Gersdorff and Borst, 2002; Zucker and Regehr, 2002). Presynaptic STD is often attributed to high pV_{r} synapses, where a single AP causes significant depletion of the RRP. In contrast, presynaptic STF has often been attributed to synapses with low initial pV_{r} and a rapid pV_{r} increase during successive APs. This was often linked to changes in Ca^{2+} signalling, for instance by rapid regulation of Ca^{2+} channels (Borst and Sakmann, 1998; Nanou and Catterall, 2018), saturation of local Ca^{2+} buffers (Eggermann et al., 2012; Felmy et al., 2003; Matveev et al., 2004), or the accumulation of intracellular Ca^{2+} which may increase pV_{r} either directly or via ‘facilitation sensors’ (Jackman and Regehr, 2017; Katz and Miledi, 1968). Alternatively, fast mechanisms increasing the RRP were proposed (Fioravante and Regehr, 2011; Gustafsson et al., 2019; Pan and Zucker, 2009; Pulido and Marty, 2017).
The coupling distance between Ca^{2+} channels and primed SVs is an important factor governing pV_{r} (Böhme et al., 2018; Eggermann et al., 2012; Stanley, 2016). Previous mathematical models describing SV fusion rates from simulated intracellular Ca^{2+} transients have in many cases relied on the assumption of uniform (or near uniform) distances between SV release sites surrounding a cluster of Ca^{2+} channels and such conditions were shown to generate STF (Böhme et al., 2016; Meinrenken et al., 2002; Nakamura et al., 2015; Vyleta and Jonas, 2014). However, alternative SV release site:Ca^{2+} channel topologies have been proposed, including two distinct perimeter distances, tight, onetoone connections of SVs and channels, or random placement of either the channels, the SVs, or both (He et al., 2019; Böhme et al., 2016; Chen et al., 2015; Guerrier and Holcman, 2018; Keller et al., 2015; Shahrezaei et al., 2006; Stanley, 2016; Wong et al., 2014). So far, the precise relationship between SV release sites and voltage gated Ca^{2+} channels on the nanometre scale is unknown for most synapses, primarily owing to technical difficulties to reliably map their precise spatial distribution. However, (M)Unc13 proteins were recently identified as a molecular marker of SV release sites (ReddyAlla et al., 2017; Sakamoto et al., 2018) and superresolution (STED) microscopy revealed that these sites surround a cluster of voltage gated Ca^{2+} channels in the center of AZs of the glutamatergic Drosophila melanogaster neuromuscular junction (NMJ) (Böhme et al., 2016; Böhme et al., 2019).
Here, by relying on the unique advantage of being able to precisely map SV release site:Ca^{2+} channel topology we study its consequence for shortterm plasticity at the Drosophila NMJ. Topologies were measured using electron microscopy (EM) following high pressure freeze fixation (HPF) or STED microscopy of Unc13 which both revealed a broad distribution of Ca^{2+} channel coupling distances. Stochastic simulations were key to identify facilitation mechanisms in the light of heterogenous SV release site:Ca^{2+} channel distances. Contrasting these simulations to physiological data revealed that models explaining STF through gradual increase in pV_{r} (from now on called ‘pV_{r}based models’) are inconsistent with the experiment while models of activitydependent regulation of the RRP account for STP profiles and synaptic variance.
Results
Distances between docked SVs and Ca^{2+} channels are broadly distributed
We first set out to quantify the SV release site:Ca^{2+} channel topology. For this we analysed EM micrographs of AZ crosssections and quantified the distance between docked SVs (i.e. SVs touching the plasma membrane) and the centre of electron dense ‘Tbars’ (where the voltage gated Ca^{2+} channels are located Fouquet et al. (2009); Kawasaki et al. (2004); Figure 1A). In wildtype animals, this leads to a broad distribution of distances (‘EM dataset wildtype’, Figure 1—figure supplement 1A; Böhme et al., 2016; Bruckner et al., 2017). At the Drosophila NMJ, the two isoforms Unc13A and –B confer SV docking and priming, but the vast majority (~95%) of neurotransmitter release and docking of SVs with short coupling distances is mediated by Unc13A (Böhme et al., 2016). We therefore investigated the docked SV distribution in flies expressing only the dominant Unc13A isoform (Unc13A rescue, see Materials and methods for exact genotypes) which showed a very similar, broad distribution of distances as wildtype animals (‘EMdataset Unc13A rescue’) (ReddyAlla et al., 2017; Figure 1A,B). In both cases, distance distributions were well described by a Rayleigh distribution (Figure 1B, Figure 1—figure supplement 1A, solid green lines). The EM micrographs studied here are a cut crosssection of a threedimensional synapse. To derive the relevant coupling distance distribution for all release sites (including the ones outside the crosssection), the Rayleigh distribution was integrated around a circle (Figure 1C), resulting in the following probability density function (pdf, see Materials and methods for derivation):
These pdfs were more symmetrical than the ones from the crosssections and peaked at larger distances (as expected from the increase in AZ area with increasing radius) (Figure 1D). The estimation of this pdf was very robust, resulting in near identical curves for the two EM datasets (Figure 1—figure supplement 1B).
We also used an independent approach to investigate the distribution of docked SV:Ca^{2+} channel coupling distances without relying on the integration of docked SV observations from crosssections: since (M)Unc13 was recently described as a molecular marker of SV release sites (ReddyAlla et al., 2017; Sakamoto et al., 2018) we investigated AZ images of wildtype NMJs stained against Unc13A (Böhme et al., 2019). Hundreds of individual AZ STED images (lateral resolution of approx. 40 nm) were aligned and averaged to obtain an average image of the AZ (Figure 1E), which revealed a ringlike distribution of the Unc13A fluorescence. In previous works we had established that the voltage gated Ca^{2+} channels reside in the center of this ring (Böhme et al., 2016). As this average image already reflects the distribution throughout the AZ area (unlike for the EM data above where an integration was necessary) the distribution of coupling distances can directly be computed based on pixel intensities and their distance to the AZ centre. Two independent datasets where analysed, resulting in very similar average images and distance distributions (‘wildtype STED dataset 1 and 2’, Figure 1—figure supplement 1).
Remarkably, although the two approaches (EM and STED microscopy) were completely independent, the distributions of coupling distances quantified by either method coincided very well (Figure 1F, Figure 1—figure supplement 1D; note that the integrated Rayleigh distributions were determined from EM micrographs and integration; they were NOT fit to the Unc13A distribution), supporting the accuracy of this realistic release site topology. The compliance between SV docking positions and Unc13A distribution further indicates that SVs dock to the plasma membrane where priming proteins are available, and therefore the entire distribution of docked SVs is potentially available for synaptic release (Imig et al., 2014).
Physiological assessment of shortterm facilitation and depression at the Drosophila NMJ
Having identified the high degree of heterogeneity in the docked SV:Ca^{2+} channel coupling distances, we became interested in how this affected synaptic function. We therefore characterized synaptic transmission at control NMJs (Ok6GAL4 crossed to w[1118]) in two electrode voltage clamp experiments. A common method to quantitatively evaluate synaptic responses and their STP behaviour is to vary the Ca^{2+} concentration of the extracellular solution which affects APinduced Ca^{2+} influx (see below). We used this approach and investigated responses evoked by repetitive (pairedpulse) AP stimulations (10 ms interval). In line with classical studies (Dodge and Rahamimoff, 1967), our results display an increase of the evoked Excitatory Junctional Current (eEJC) responses to the first AP (eEJC_{1} amplitudes) with increasing extracellular Ca^{2+} (Figure 2A,B). STP was assessed by determining the pairedpulse ratio (PPR): the amplitude of the second response divided by first. The eEJC_{2}amplitude was determined taking the decay of eEJC_{1} into account (see insert in Figure 2C, Figure 2—figure supplement 1A). At low extracellular Ca^{2+} (0.75 mM), we observed strong STF (with an average PPR value of 1.80), which shifted towards depression (PPR < 1) with increasing Ca^{2+} concentrations (Figure 2C,D). Thus, the same NMJ displays both facilitation and depression depending on the extracellular Ca^{2+} concentration, making this a suitable model synapse to investigate STP behaviour.

Figure 2—source data 1
 https://cdn.elifesciences.org/articles/51032/elife51032fig2data1v2.zip

Figure 2—source data 2
 https://cdn.elifesciences.org/articles/51032/elife51032fig2data2v2.zip

Figure 2—source data 3
 https://cdn.elifesciences.org/articles/51032/elife51032fig2data3v2.xlsx
In panels B and D the mean eEJC_{1} amplitudes and PPRs from six animals are shown and the error bars indicate standard deviation, SD (across all animals). We also examined the variation of repeated APevoked responses at the same NMJ between trials (10 s apart) at different extracellular Ca^{2+} concentrations (Figure 2E,F). At low concentrations (0.75 mM), the probability of transmitter release is low, resulting in a low mean eEJC_{1} amplitude with little variation (Figure 2E,F, Figure 2—figure supplement 2 ). With increasing extracellular Ca^{2+}, the likelihood of SV fusion increased and initially so did the variance (e.g. at 1.5 mM extracellular Ca^{2+}). However, further increase in extracellular Ca^{2+} (3 mM, 6 mM, 10 mM) led to a drop in variance (Figure 2E, Figure 2—figure supplement 2). Figure 2F depicts this average ‘variancemean’ relationship from 6 cells (means of cell means and means of cell variances, error bars indicate SEM). When assuming a binomial model, this approach has often been used to estimate the number of release sites n_{sites} and the size of the postsynaptic response elicited by a single SV (q) (Clements and Silver, 2000). In agreement with previous studies of the NMJ this relationship was well described by a parabola with forced intercept at y = 0 and n_{sites} = 164 and q = 0.64 nA (Figure 2F, Figure 2—figure supplement 2; Matkovic et al., 2013; Müller et al., 2012; Weyhersmüller et al., 2011).
Simulation of APinduced Ca^{2+} signals
Having determined the distribution of coupling distances (Figure 1) and the physiological properties of the NMJ synapse (Figure 2), we next sought to compare how the one affected the other. There are two things two consider here. First of all, the SV release probability steeply depends on the 4^{th} to 5^{th} power of the local Ca^{2+} concentration (Neher and Sakaba, 2008). Secondly, because of the strong buffering of Ca^{2+} signals at the synapse, the magnitude of the APevoked Ca^{2+} transients dramatically declines with distance from the Ca^{2+} channel (Böhme et al., 2018; Eggermann et al., 2012). These two phenomena together make the vesicular release probability extremely sensitive to the coupling distance to the Ca^{2+} channels. Because we find that this distance is highly heterogeneous among SVs within the same NMJ, the question arises how these two properties (heterogeneity of distances combined with a strong distance dependence of pV_{r}) functionally impact on synaptic transmission. Indeed, approaches by several labs to map the activity of individual NMJ AZs revealed highly heterogeneous activity profiles (Akbergenova et al., 2018; Gratz et al., 2019; Muhammad et al., 2015; Peled and Isacoff, 2011).
To quantitatively investigate the functional impact of heterogeneous SV placement, we wanted to use mathematical modelling to predict APinduced fusion events of docked SVs placed according to the found distribution. A prerequisite for this is to first faithfully simulate local, APinduced Ca^{2+} signals throughout the AZ (such that the local transients at each docking site are known). We first determined the relevant AZ dimensions at the Drosophila NMJ, which, similarly to the murine Calyx of Held, is characterized by many AZs operating in parallel. We therefore followed previous suggestions from the Calyx using a box with reflective boundaries containing a cluster of Ca^{2+} channels in the base centre (Meinrenken et al., 2002). The base dimensions (length = width) were determined as the mean interAZ distance of all AZs to their four closest neighbours (because of the 4fold symmetry) from NMJs stained against the AZmarker BRP (Kittel et al., 2006; Wagh et al., 2006; Figure 3A). To save computation time, we further simplified to a cylindrical simulation (where the distance to the Ca^{2+} channel is the only relevant parameter) covering the same AZ area (Figure 3B, Table 1).
To simulate the electrophysiological experiments above, where the extracellular Ca^{2+} concentration was varied (Figure 2), it was important to establish how the extracellular Ca^{2+} concentration influenced APinduced Ca^{2+} influx. In particular, it is known that Ca^{2+} currents saturate at high extracellular Ca^{2+} concentrations (Church and Stanley, 1996). Unlike other systems, the presynaptic NMJ terminals are not accessible to electrophysiological recordings, so we could not measure the currents directly. We therefore used a fluorescencebased approach as a proxy. APevoked Ca^{2+} influx was assessed in flies presynaptically expressing the Ca^{2+}dependent fluorescence reporter GCaMP6m (;P{y[+t7.7] w[+mC]=20XUASIVSGCaMP6m}attP40/Ok6GAL4). Fluorescence increase was monitored upon stimulation with 20 APs (at 20 Hz) while varying the extracellular Ca^{2+} concentration and showed saturation behaviour for high concentrations (Figure 3—figure supplement 1). This is consistent with a previously described MichaelisMenten type saturation of fluorescence responses of a Ca^{2+}sensitive dye upon single AP stimulation at varying extracellular Ca^{2+} concentrations at the Calyx of Held, where halfmaximal Ca^{2+} influx was observed at 2.6 mM extracellular Ca^{2+} (Schneggenburger et al., 1999). This relationship was successfully used in the past to predict Ca^{2+} influx in modeling approaches Trommershäuser et al. (2003). In our measurements, we determined a half maximal fluorescence response at a very similar concentration of 2.68 mM extracellular Ca^{2+} and therefore used this value as K_{M,current} in a MichaelisMenten equation (Materials and methods, Equation 5) to calculate APinduced presynaptic Ca^{2+} influx. The second parameter of the MichaelisMenten equation, (the maximal Ca^{2+} current charge, Q_{max}) was optimized for each model (Figure 3—figure supplement 2, for parameter explanations and best fit parameters see Table 2). We furthermore assumed that basal, intracellular Ca^{2+} concentrations at rest were also slightly dependent on the extracellular Ca^{2+} levels in a MichaelisMenten relationship with the same dependency (K_{M,current}) and a maximal resting Ca^{2+} concentration of 190 nM (resulting in 68 nM presynaptic basal Ca^{2+} concentration at 1.5 mM external Ca^{2+}). With these and further parameters taken from the literature on Ca^{2+} diffusion and buffering (see Table 1) the temporal profile of Ca^{2+} signals in response to paired AP stimulation (10 ms interval) could be calculated at all AZ locations using the software CalC (Matveev et al., 2002; Figure 3C, Figure 3—figure supplement 2). This enabled us to perform simulations of NT release from vesicles placed according to the distribution described above.
Stochastic simulations and fitting of release models
In the past, we and others have often relied on deterministic simulations based on numerical integration of kinetic reaction schemes (ordinary differential equations, ODEs). These are computationally effective and fully reproducible, making them wellbehaved and ideal for the optimisation of parameters (a property that was also used here for initial parameter searches, see Materials and methods). However, NT release is quantal and relies on only a few (hundred) SVs, indicating that stochasticity plays a large role (Gillespie, 2007). Moreover, deterministic simulations always predict identical output making it impossible to analyse the synaptic variance between successive stimulations, which is a fundamental hallmark of synaptic transmission and an important physiological parameter (Figure 2F; Scheuss and Neher, 2001; VereJones, 1966; Zucker, 1973). Stochastic simulations allow a prediction of variance which can help identify adequate models that will not only capture the mean of the data, but also its variance. To compare this, data points are now shown with error bars indicating the square root of the average variance between stimulations within a cell (Figure 4C, E, 6E, G and 7E, G). This is the relevant parameter since the model is designed to resemble an ‘average’ NMJ’ and therefore cannot predict interanimal variance. Finally, as we show here deterministic simulations cannot be compared to experimentally determined PPR values because of Jensen’s inequality (full proof in Materials and methods, see Figure 4—figure supplement 1). Thus, stochastic simulations are necessary to account for SV pool sizes, realistic release site distributions, synaptic variance and STP. We thus implemented stochastic models of SV positions (drawn randomly from the distribution above) and SV Ca^{2+} binding states based on inhomogeneous, continuous time Markov models with transition rates governing reaction probabilities (see Materials and methods for details).
We also needed to consider where new SVs would (re)dock once SVs had fused and implemented the simplest scenario of redocking in the same positions. This ensures a stable distribution over time and agrees with the notion that vesicles prime into predefined release sites, which are stable over much longer time than a single priming/unpriming event (ReddyAlla et al., 2017).
A singlesensor model fails to induce sufficient facilitation and produces excessive variance
The first model we tested was the singlesensor model proposed by Lou et al. (2005), where an SV binds up to 5 Ca^{2+} ions, with each ion increasing its fusion rate or probability (Figure 4A, Table 3). Release sites were placed according to the distance distribution in Figure 1D and all sites were occupied by a primed SV prior to stimulation (i.e. the number of release sites equals the number of vesicles in the RRP). Sites becoming available following SV fusion were replenished from an unlimited vesicle pool, making the model identical to the one described by Wölfel et al. (2007). Ca^{2+} (un)binding kinetics were taken from Wölfel et al. (2007) Table 3, the values of the maximal Ca^{2+} current charge (Q_{max}), the SV replenishment rate (k_{rep}) and the number of release sites (n_{sites}) were free parameters optimized to match the experimental data (see Materials and methods for details, best fit parameters in Table 2).
To be able to compare the output of this and all subsequent models to experimental data as depicted in Figure 2 (postsynaptic eEJC measurements), the predicted fusion events were convolved with a typical postsynaptic response to the fusion of a single SV (mEJC, Figure 2—figure supplement 1B, see Materials and methods for more details). From the stochastic simulations (1000 runs each), we calculated the mean and variance of eEJC_{1} amplitudes, and the mean and variance of PPRs at various extracellular Ca^{2+} concentrations and contrasted these with the experimental data.
This singlesensor model was able to reproduce the eEJC_{1} values (Figure 4B,C). Moreover, the model accounted for the STD typically observed at high extracellular Ca^{2+} concentrations in the presence of rapid replenishment (Hallermann et al., 2010; Miki et al., 2016) (our best fit yielded τ ≈ 6 ms and reducing this rate led to unnaturally strong depression, Figure 4E, green curve+area). However, even despite rapid replenishment this model failed to reproduce the STF observed at low extracellular Ca^{2+} (Figure 4D,E) and the variances predicted by this model were much larger than found experimentally (Figure 4F,G). The observation that eEJC_{1} amplitudes were well accounted for, but STPs were not, may relate to the fact that this model was originally constructed to account for a single Ca^{2+}triggered release event (Lou et al., 2005). As this model lacks a specialized mechanism to induce facilitation, residual Ca^{2+} binding to the Ca^{2+ }sensor is the only facilitation method which appears to be insufficient (Jackman and Regehr, 2017; Ma et al., 2015; Matveev et al., 2002). This result differs from our previous study using this model where we had placed all SVs at the same distance to Ca^{2+} channels which reliably produced STF (Böhme et al., 2016). So why does the same model fail to produce STF with this broad distribution of distances? To understand this we investigated the spatial distribution of SV release in simulations of the pairedpulse experiment at 0.75 mM extracellular Ca^{2+} (Figure 5).
In the absence of a facilitation mechanism, only part of the SV distribution is utilized
Figure 5A depicts two examples of synapses – seen from above – with SVs randomly placed according to the distance distribution in Figure 1D/5B. The synapse is shown immediately before AP_{1}, immediately after AP_{1}, immediately before AP_{2} (i.e. after refilling) and immediately after AP_{2} (the external Ca^{2+} concentration was 0.75 mM). From this analysis it becomes clear that the pV_{r}1 caused by AP_{1} essentially falls to zero around the middle of the SV distribution (Figure 5B, top panel). This means that only SVs close to the synapse center fuse, and these highpV_{r} SVs are depleted by AP_{1}. SV replenishment refills the majority (but not all) of those sites and thus AP_{2}/pV_{r}2 essentially draws on the same part of the distribution (Figure 5B, bottom panel). Because of this, and because the refilling is incomplete, this causes STD. Even with faster replenishment (which would be incompatible with the low PPR values at high extracellular Ca^{2+}, Figure 4E) this scenario would only lead to a modest increase of the PPR to values around 1. Therefore, our analysis reveals that large variation in Ca^{2+} channel distances results in a specific problem to generate STF. Our analysis further indicates that with the best fit parameters of the single sensor model, the majority of SVs (those further away) is not utilized at all.
A dual fusionsensor model improves PPR values, but generates too little facilitation and suffers from asynchronous release and too much variance
The singlesensor model failed to reproduce the experimentally observed STF at low extracellular Ca^{2+} concentrations because of the dominating depletion of SVs close to Ca^{2+} channels, and the inability to draw on SVs further away. However, this situation may be improved by a second Ca^{2+} sensor optimized to enhance the pV_{r}2 in response to AP_{2}. Indeed, in the absence of the primary Ca^{2+} sensor for fusion, Ca^{2+} sensitivity of synaptic transmission persists, which was explained by a dual sensor model (Sun et al., 2007). It was recently suggested that syt7 functions alongside syt1 as a Ca^{2+} sensor for release (Jackman et al., 2016), and deterministic mathematical dual fusionsensor model assuming homogeneous release probabilities (which implies homogeneous SV release site:Ca^{2+} channel distances) was shown to generate facilitation (Jackman and Regehr, 2017). Similarly, stochastic modelling of NT release at the frog NMJ also showed a beneficial effect of a second fusion sensor for STF (Ma et al., 2015). We therefore explored whether a dual fusion sensor model could account for synaptic facilitation from realistic release site topologies.
The central idea of this dual fusionsensor model is that while syt1 is optimized to detect the rapid, APinduced Ca^{2+} transients (because of its fast Ca^{2+} (un)binding rates, but fairly low Ca^{2+} affinity), the cooperating Ca^{2+} sensor is optimized to sense the residual Ca^{2+} after this rapid transient (Figure 3C) (with slow Ca^{2+} (un)binding, but high Ca^{2+ }affinity). The activation of this second sensor after (but not during) AP_{1} could then enhance the release probability of the remaining SVs for AP_{2} (Figure 6A,B). This is illustrated in Figure 6B, where k_{2} (the onrate of Ca^{2+} binding to the slow sensor) is varied resulting in different time courses and amounts of Ca^{2+} binding to the second sensor. Increasing the release probability is equivalent to lowering the energy barrier for SV fusion (Schotten et al., 2015). In this model both sensors regulate pV_{r} and therefore additively lower the fusion barrier with each associated Ca^{2+} ion (Figure 6A), resulting in multiplicative effects on the SV fusion rate. While the fast fusion reaction appears to have a 5fold Ca^{2+} cooperativity (Bollmann et al., 2000; Burgalossi et al., 2010; Schneggenburger and Neher, 2000), it is less clear what the Ca^{2+} cooperativity of a second Ca^{2+} sensor may be, although the fact that the cooperativity is reduced in the absence of the fast sensor (Burgalossi et al., 2010; Kochubey and Schneggenburger, 2011; Sun et al., 2007) could be taken as evidence for a Ca^{2+} cooperativity < 5. We explored cooperativities 2, 3, 4, and 5 (cooperativities 2 and 5 are displayed in Figure 6 and Figure 6—figure supplement 1). It is furthermore not clear whether such a sensor would be targeted to the SV (like syt1 /2), or whether it is present at the plasma membrane. Both scenarios are functionally possible and it was indeed reported that syt7 is predominantly or partly localized to the plasma membrane (Sugita et al., 2001; Weber et al., 2014). A facilitation sensor on the plasma membrane would be more effective, which our simulations confirmed (not shown), because it would not be consumed by SV fusion, allowing the sensor to remain activated. We therefore present this version of the model here. We used a second sensor with a Ca^{2+} affinity of K_{D} = 1.5 μM (Brandt et al., 2012; Jackman and Regehr, 2017).
Like for the singlesensor model, all release sites were occupied with releasable vesicles (n_{sites} equals the number of RRP vesicles) and their locations determined by drawing random numbers from the pdf. When fitting this model five parameters were varied: Q_{max}, k_{rep}, and n_{sites} (like in the singlesensor model) together with k_{2} (Ca^{2+} association rate constant to the second sensor) and s (the factor describing the effect of the slow sensor on the energy barrier for fusion) (see Table 2 for best fit parameters). The choice of k_{2} had an effect on the PPR in simulations, confirming that the second sensor was able to improve the release following AP_{2} (Figure 6C). Figure 6D–I show that the dual fusionsensor model could fit the eEJC_{1} amplitudes and the model slightly improved the higher PPR values at the low and the lower PPR values at high extracellular Ca^{2+} concentrations compared to the single sensor model (compare Figures 4E and 6G). However, the model failed to produce the STF observed experimentally (the PPR values at 0.75 mM Ca^{2+}were ~ 1.08 in the simulation compared to ~ 1.80 in the experiments). Another problem of the dual fusionsensor model was that release became more asynchronous than observed experimentally (Figure 6D), which was due to the triggering of SV fusion inbetween APs. Finally, predicted variances were much larger than the experimental values (Figure 6I).
In addition to the optimization, we systematically investigated a large region of the parameter space (Figure 6J,K), but found no combination of parameters that would be able to generate the experimentally observed STF. Lowering the Ca^{2+} influx (by decreasing Q_{max}) yielded a modest increase in PPR values (Figure 6J), but required a large number of release sites (n_{sites}) to match the eEJC_{1} amplitudes (Figure 6K). Changing s had the largest effect when k_{2} was close to the best fit value and moving away from this value decreased the PPRs, either by increasing the effect of the second sensor on AP_{1} (when increasing k_{2}) or by decreasing the effect on AP_{2} (when decreasing k_{2}), which both counteracts STF (Figure 6B,J).
Fitting the dual fusionsensor model with a Ca^{2+} cooperativity of 5 did not improve the situation (Figure 6—figure supplement 1, best fit parameters in Table 2): Although slightly more facilitation was observed, this model suffered from even larger variance overshoots (Figure 6—figure supplement 1E) and excessive asynchronous release (Figure 6—figure supplement 1A,C). We explored different K_{D} values between 0.5 and 2 μM at cooperativities 2–5 in separate optimizations, but found no satisfactory fit of the data (results not shown). Thus, a dual fusionsensor model is unlikely to account for STF observed from the realistic SV release site topology at the Drosophila NMJ. Note that this finding does not rule out that syt7 functions in STF, but argues against a role in cooperating alongside syt1 in a pV_{r}based facilitation mechanism.
Rapidly regulating the number of RRP vesicles accounts for eEJC_{1} amplitudes, STF, temporal transmission profiles and variances
Since dual fusionsensor models and other models depending on changes in pV_{r} (see Discussion) are unlikely to be sufficient, we next investigated mechanisms involving an activitydependent regulation of the number of participating release sites. For this we extended the singlesensor model by a single unpriming reaction (compare Figures 4A and 7A). The consequence of reversible priming is that the initial release site occupation can be less than 100% (in which cases n_{sites} can exceed the number of RRP vesicles). This enables an increase (‘overfilling’) of the RRP (/increase in site occupancy) during the interstimulus interval (consistent with reports in other systems Dinkelacker et al., 2000; Gustafsson et al., 2019; Pulido et al., 2015; Smith et al., 1998; Trigo et al., 2012). We assumed that Ca^{2+} would stabilize the RRP/release site occupation by slowing down unpriming (Figure 7A). This made the steadystate RRP size dependent on the resting Ca^{2+} concentration and the modest dependence of this on the extracellular Ca^{2+} resulted in RRP enlargement with increasing extracellular Ca^{2+} (Figure 7B), in agreement with recent findings on central synapses (Malagon et al., 2020). This model (like the dual fusionsensor models depicted in Figure 6 and Figure 6—figure supplement 1) includes two different Ca^{2+ }sensors, but the major difference is that these Ca^{2+} sensors operate to regulate two separate sequential steps (priming and fusion). Indeed, this scenario aligns with reports of a syt7 function upstream of SV fusion (Liu et al., 2014; Schonn et al., 2008). Figure 7C shows how the number of RRP vesicles develops over time in this model during a pairedpulse experiment for low and high extracellular Ca^{2+} concentrations. In all cases, SV priming was in equilibrium prior to the first stimulus, indicated by the horizontal lines (0–2 ms, Figure 7C). Note that prior to AP_{1} priming is submaximal (~41%) for 0.75 mM extracellular Ca^{2+}, but near complete (~99%) at 10 mM extracellular Ca^{2+}. At low extracellular Ca^{2+} the elevation of Ca^{2+} caused by AP_{1} results in a sizable inhibition of unpriming, leading to an increase (‘overfilling’) of the RRP during the interstimulus interval. With this, more primed SVs are available for AP_{2}, causing facilitation (green line in Figure 7C). In contrast, at high extracellular Ca^{2+} concentrations, the rate of unpriming is already low at steady state and the RRP close to maximal capacity (grey line in Figure 7C). At this high extracellular Ca^{2+} concentration, AP_{1} induces a larger Ca^{2+} current (higher pV_{r}), resulting in strong RRP depletion, of which only a fraction recovers between APs (as in the other models, replenishment commences with a Ca^{2+} independent rate k_{rep}). Because Ca^{2+} acts in RRP stabilization, not in stimulating forward priming, this model (unlike the dual fusionsensor models in Figure 6 and Figure 6—figure supplement 1) did not yield asynchronous release inbetween APs (Figure 7D). Thus, the two most important features of this model are the submaximal site occupation and an inhibition of unpriming by intracellular Ca^{2+}.
In this model we assumed a Ca^{2+} cooperativity of n = 5 for the unpriming mechanism (we also explored n = 2, see Figure 7—figure supplement 1). The following parameters were optimized: Q_{max}, n_{sites} and k_{rep} (like in the single and dual fusionsensor models), together with K_{M,prim}, the Ca^{2+} affinity of the priming sensor, and u, its Ca^{2+} cooperativity. These values together define the Ca^{2+}dependent unpriming rate (see Table 2 for best fit parameters). The total number of fitted parameters (5) was the same as for the dual fusionsensor models (Figure 6 and Figure 6—figure supplement 1). Figure 7D–I present the results. It is clear that both eEJC_{1} amplitudes and PPR values were described very well with this model at all extracellular Ca^{2+} concentrations. In addition, the variancemean relationship of the eEJC_{1} was reproduced satisfactorily, except for a small variance overshoot for the highest extracellular Ca^{2+} concentrations (Figure 7I, see Discussion). Fitting of the unpriming model with a Ca^{2+} cooperativity of 2 also led to a good fit (Figure 7—figure supplement 1), although the variance overshoot was somewhat larger. We also explored the timedependence of the facilitation by simulating PPR values for various interstimulus intervals at different extracellular Ca^{2+} concentrations which could be investigated experimentally in the future to further refine parameters (Figure 7—figure supplement 2).
Different facilitating synapses exhibit a large range of PPR values, some larger than observed at the Drosophila NMJ (Jackman et al., 2016). Therefore, if this were a general mechanism to produce facilitation, we would expect it to be flexible enough to increase the PPR much more than observed here. To investigate the model’s flexibility we systematically explored the parameter space by varying Q_{max}, K_{M,prim}, and u (Figure 7J,K). Similar to Figure 6J,K, the colors of the balls represent the PPR value and the number of release sites needed to fit the eEJC_{1} amplitudes. Consistent with a very large dynamic range of this mechanism, PPR values ranged from 0.85 to 3.90 (Figure 7J,K) and unlike the dual fusionsensor model, PPR values were fairly robust to changes in Ca^{2+} influx (note the different scales on Figure 7J,K and Figure 6J,K). Moreover, because this mechanism does not affect the Ca^{2+} sensitivity of SV fusion, facilitation was achieved without inducing asynchronous release (Figure 7D).
We also investigated an alternative model based on Ca^{2+}dependent release site activation. In this model, all sites are occupied by a vesicle, but some sites are inactive and fusion is only possible from activated sites. We assumed that site activation was Ca^{2+}dependent. In order to avoid site activation during AP_{1}, which would again hinder STF and could contribute to asynchronous release, we implemented an intermediate delay state (Figure 7—figure supplement 3A–B) from which sites were activated in a Ca^{2+}independent reaction. This could mean that priming occurs in twosteps, with the first step being Ca^{2+}dependent. Similar to the unpriming model presented above, the modest increase of intracellular Ca^{2+} with extracellular Ca^{2+} yielded an RRP increase (/increase in active sites) (Figure 7—figure supplement 3I). This model agreed similarly well with the data as the unpriming model (Figure 7—figure supplement 3C–H). Thus, both mechanisms which modulate the RRP rather than pV_{r} are fully capable of reproducing the experimentally observed Ca^{2+}dependent eEJC_{1} amplitudes, STF, release synchrony and variance. The unpriming model was preferred since it had fewer parameters and performed slightly better in optimisations than the site activation model.
A release site facilitation mechanism utilizes a larger part of the SV distribution
Why do n_{site}/primingbased mechanisms (Figure 7, Figure 7—figure supplement 1, Figure 7—figure supplement 3) account for STF from the broad distribution of SV release site:Ca^{2+} channel coupling distances, while the pV_{r}based models (Figures 4 and 6, Figure 6—figure supplement 1) cannot? To gain insight into this, we analysed the spatial dependence of transmitter release in the unpriming model during the pairedpulse experiment (0.75 mM extracellular Ca^{2+}) in greater detail (Figure 8). Panel 8A, similarly to Figure 5A, shows example stochastic simulations (at external Ca^{2+} concentration 0.75 mM, to illustrate facilitation). The best fit parameters of the unpriming model predicted a larger Ca^{2+} influx (1.64fold and 3.05fold larger Q_{max} value) than the single and dual fusionsensor models (Table 2). The larger Ca^{2+} influx compensated for the submaximal priming of SVs (reduced release site occupancy) prior to the first stimulus by expanding the region where SVs are fused (Figure 8B). Comparing to Figure 5B, a much larger part of the SV distribution is utilized during the first stimulus. Following AP_{1}, vesicles prime into empty sites across the entire distribution, allowing AP_{2} to draw again from the entire distribution. During this time, the increased residual Ca^{2+} causes overfilling of the RRP, that is more release sites are now occupied, giving rise to more release during AP_{2}. Notably, the AP_{2}induced release again draws from the entire distribution. Thus, the unpriming model not only reproduces STF and synaptic variance, but also utilizes docked SVs more efficiently from the entire distribution compared to the single and dual fusionsensor model.
Discussion
We here described a broad distribution of SV release site:Ca^{2+} channel coupling distances in the Drosophila NMJ and compared physiological measurements with stochastic simulations of four different release models (singlesensor, dual fusionsensor, Ca^{2+}dependent unpriming and site activation model). We showed that the two first models (singlesensor and dual fusionsensor), where residual Ca^{2+} acts on the energy barrier for fusion and results in an increase in pV_{r}, failed to reproduce facilitation. The two latter models involve a Ca^{2+}dependent regulation of participating release sites and reproduced release amplitudes, variances and PPRs. Therefore, the Ca^{2+}dependent accumulation of releasable SVs is a plausible mechanism for pairedpulse facilitation at the Drosophila NMJ, and possibly in central synapses as well. In more detail, our insights are as follows:
The SV distribution was described by the singlepeaked integrated Rayleigh distribution with a fitted mean of 122 nm. The distribution has a low probability for positioning of SVs very close to Ca^{2+ }channels (less than 1.5% within 30 nm) and is therefore reasonably consistent with suggestions of a SV exclusion zone of ~ 30 nm around Ca^{2+ }channels (Keller et al., 2015). Strikingly, almost exactly the same distribution was identified for the essential priming protein Unc13A (Figure 1F, Figure 1—figure supplement 1D), indicating that docked SVs are likely primed (Imig et al., 2014).
The broad distribution of SV release site:Ca^{2+} channel distances particularly impedes pV_{r}based facilitation mechanisms. Indeed, previous models that reproduced facilitation using pV_{r}mechanisms typically placed SVs at an identical/similar distance to Ca^{2+} channels, resulting in intermediate (and identical) pV_{r} for all SVs (Böhme et al., 2016; Böhme et al., 2018; Bollmann and Sakmann, 2005; Fogelson and Zucker, 1985; Jackman and Regehr, 2017; Matveev et al., 2006; Matveev et al., 2004; Tang et al., 2000; Vyleta and Jonas, 2014; Yamada and Zucker, 1992). Here, having mapped the precise AZ topology, we show that the broad SV distribution together with the steep dependence of release rates on [Ca^{2+}] creates a situation where pV_{r} falls to almost zero for SVs further away than the mean of the distribution (Figure 5). As a result, most SVs either fuse during AP_{1}, or have pV_{r} values close to zero, leaving little room for modulation of pV_{r} to create facilitation. Such mechanisms (including buffer saturation, and Ca^{2+ }binding to a second fusion sensor) will act to multiply release rates with a number > 1. However, since SVs with pV_{r} close to one have already fused during AP_{1}, and most of the remaining vesicles have pV_{r} close to zero such a mechanism will be ineffective in creating facilitation. Thus, the broad distribution of SV release site:Ca^{2+} channel distances makes it unlikely that pV_{r}based mechanisms can cause facilitation.
The dual fusionsensor model was explored as an example of a pV_{r}based model. Two problems were encountered: The first problem was that the second sensor, due to its high affinity for Ca^{2+}, was partly activated in the steady state prior to the stimulus (Figure 6B). Therefore, it could not increase pV_{r}2 without also increasing pV_{r}1. This makes it inefficient in boosting the PPR. The second problem was kinetic: the second sensor should be fast enough to activate between two APs, but slow enough not to activate during AP_{1}. This is illustrated in Figure 6B–C, which shows the time course of activation of the two sensors and the corresponding PPR values for varying Ca^{2+} binding rates of the second sensor. Since the sensor is Ca^{2+}dependent, the rate inevitably increases during the Ca^{2+} transient, leading to too much asynchronous release. In principle, the first problem could be alleviated by increasing the Ca^{2+ }cooperativity of the second sensor, which would make it easier to find parameters where the sensor would activate after but not before AP_{1}. We therefore tried to optimize the model with cooperativities of 3, 4, and 5 (Figure 6—figure supplement 1 shows cooperativity 5), and indeed, the higher cooperativity made it possible to obtain slightly more facilitation. However, activation during the AP (the second problem) was exacerbated and caused massive and unphysiological asynchronous release. Thus, a secondary Ca^{2+} sensor acting on the energy barrier for fusion is unlikely to account for facilitation in synapses with a broad distribution of SV release site:Ca^{2+} channel distances.
We included stochasticity at the level of the SV distribution (release sites were randomly drawn from the distribution) and at the level of SV Ca^{2+} (un)binding and fusion. This was essential since deterministic and stochastic simulations do not agree on PPRvalues due to Jensen’s inequality (for a stochastic process the mean of a ratio is not the same as the ratio of the means) (see Materials and methods and Figure 4—figure supplement 1). The effect is largest when the evoked release amplitude is smallest. Since small amplitudes are often associated with high facilitation, this effect is important and needs to be taken into account. Stochastic Ca^{2+ }channel gating on the other hand was not included, as this would increase simulation time dramatically. At the NMJ, the Ca^{2+} channels are clustered (Gratz et al., 2019; Kawasaki et al., 2004), and most SVs are relatively far away from the cluster, a situation that was described to make the contribution of Ca^{2+ }channel gating to stochasticity small (Meinrenken et al., 2002). However, the situation will be different in synapses where individual SVs colocalize with individual Ca^{2+ }channels (Stanley, 2016).
Stochastic simulations made it possible to not only determine the mean eEJC_{1} and PPR values, but also the standard deviation around these values upon repeated activation of the NMJ (indicated as lightly colored bands on the simulations in Figure 4C,E, Figure 6E,G, and Figure 7E,G), which can be compared to measurements (shown as black error bars in the same figure panels). This also enabled us to compare our model to experimental variancemean data (Figures 4G, 6I and 7I), which we found was key to identify valid models. All models tested resulted in variancemean dependences that were well approximated by a parabola with intercept 0. Note that such parabola agrees with the meanvariance relationship in a binomial distribution. However these simulations show that the assumption of heterogeneous release probability (and changing RRP size) can also lead to the experimentally observed parabolic variancemean relationship. The singlesensor and dual fusionsensor models resulted in overshooting variances (Figures 4G and 6I), which became even worse in the case of higher cooperativity of the second fusion sensor (Figure 6—figure supplement 1). The righthand intercept of the variancemean relationship with the abscissa is interpreted as the product of the number of release sites (n_{sites}) and q (the single SV quantum) and the tendency of these models to overshoot the variance is due to the fitting procedure increasing n_{sites}, while at the same time reducing pV_{r}, (by reducing Q_{max}, the maximal APinduced Ca^{2+ }influx). The lower pV_{r} increases the PPR by reducing the effect of depletion, but results in unrealistically high n_{sites.} Therefore, it was essential to contrast the models to experimental variancemean data, which restrict n_{sites}. This revealed that pV_{r}based facilitation mechanisms produced unrealistic variancemean behavior. In this context, models involving a Ca^{2+}dependent accumulation of releasable SVs fare much better, because only those can cause facilitation in the presence of realistic n_{sites}, resulting in very similar variancemean behaviour to the experiment (Figure 7I, Figure 7—figure supplement 1E). The remaining slight overshoot for variances at high extracellular Ca^{2+} concentrations could have technical/experimental reasons, because these experiments are of long duration, which might lead to rundown over time (which is not present in the model simulations) that causes a compression of the parabolic relationship along the abscissa (experiments were performed by increasing Ca^{2+} concentrations).
We arrived at two models that can explain pairedpulse facilitation and variancemean behaviour at the Drosophila NMJ. Both models include a Ca^{2+}dependent increase in the number of participating (occupied/activated) release sites. In the Ca^{2+}dependent unpriming model, forward priming happens at a constant rate, but unpriming is inversely Ca^{2+}dependent, such that increases in residual Ca^{2+} lead to inhibition of unpriming, thereby increasing release site occupation between stimuli (Figure 7). Ca^{2+}dependent replenishment has been observed in multiple systems (Dinkelacker et al., 2000; Smith et al., 1998; Stevens and Wesseling, 1998; Wang and Kaczmarek, 1998). This has traditionally been implemented in various release models as a Ca^{2+}dependent forward priming rate (Man et al., 2015; Pan and Zucker, 2009; Voets, 2000; Weis et al., 1999). In a previous secretion model in chromaffin cells, we had proposed a catalytic function of Ca^{2+} upstream of vesicle fusion (Walter et al., 2013). However, in the context of STF such models would favour accelerated priming during the AP, which would counteract this facilitation mechanism and might cause asynchronous release, similar to the problem with the dual fusionsensor model (Figure 6). In the model presented here this is prevented by including the Ca^{2+} dependency on the unpriming rate. Consistent with this idea, recent data in cells and in biochemical experiments showed that the Ca^{2+}dependent priming protein (M)Unc13 reduces unpriming (He et al., 2017; Prinslow et al., 2019). Another model that reproduced the electrophysiological data was the site activation model, where sites are activated Ca^{2+}dependently under docked (but initially unprimed) SVs (Figure 7—figure supplement 3). In this case, we had to prevent rapid activationandfusion during the AP by including an extra, Ca^{2+}independent transition, which introduces a delay before sites are activated (Figure 7—figure supplement 3). The two models are conceptually similar in that they either recruit new SVs to (always active) sites, or activate sites underneath dormant SVs. Those two possibilities are almost equivalent when measuring with electrophysiology, but they might be distinguished in the future using flashandfreeze electron microscopy (Chang et al., 2018; Watanabe et al., 2013). Interestingly, Unc13 has recently been shown to form release sites at the Drosophila NMJ (ReddyAlla et al., 2017). Therefore, the two models also correspond to two alternative interpretations of Unc13 action (to prevent unpriming, or form release sites).
In our model, all primed vesicles have identical properties, and only deviate in their distance to the Ca^{2+} channel cluster (positional priming, Neher and Brose, 2018). Alternatively, several vesicle pools with different properties (molecular priming) could be considered, which might involve either vesicles with alternative priming machineries, or vesicles being in different transient states along the same (slow) priming pathway (Walter et al., 2013). In principle, if different primed SV states are distributed heterogenously such that more distant vesicles are more primed/releasable, such an arrangement might counteract the effects of a broad distance distribution, although this is speculative. Without such a peripheral distribution, the existence of vesicles in a highly primed/releasable state (such as the ‘superprimed’ vesicles reported at the Calyx of Held synapse), would result in pronounced STD, and counteract STF, which indeed has been observed (Lee et al., 2013; Taschenberger et al., 2016).
In this study electrophysiological recordings were performed on muscle 6 of the Drosophila larva which receives input from morphologically distinct NMJs containing big (Ib) and small (Is) synaptic boutons, which have been shown to differ in their physiological properties (Atwood et al., 1993; He et al., 2009; Newman et al., 2017). This could add another layer of functional heterogeneity in the postsynaptic responses analysed here (the EM and STED analyses shown here were focused on Ib inputs). Because our model does not distinguish between Is and Ib inputs, the estimated parameters represent a compound behaviour of all types of synaptic input to this muscle. Future investigations to isolate the contribution of the different input types (e.g. by genetically targeting Is/Ibspecific motoneurons using recently described GAL4 lines; PérezMoreno and O'Kane, 2019) could help distinguish between inputs and possibly further refine the model to identify parameter differences between these input types.
Figure 9 summarizes the results for the singlesensor, dual fusionsensor and unpriming models. Facilitation in single and dual fusionsensor models depend on the increase in release probability from the first AP to the next (compare colored rings representing 25% release probability between row 2 and 4). However, the increase is very small, even for the dual fusionsensor model, and to nevertheless produce some facilitation, optimisation finds a small Ca^{2+ }influx, which leads to an ineffective use of the broad vesicle distribution (and a toohigh estimate of n_{sites}). In the unpriming model a higher fitted Ca^{2+} influx (Q_{Max}) leads to a more effective use of the entire SV distribution, and facilitation results from the combination of incomplete occupancy of release sites before the first AP (row 1), combined with ‘overshooting’ priming into empty sites between APs (row 3).
Molecularly, syt7 was linked to STF behaviour (Jackman et al., 2016), and our data does not rule out that syt7 is essential for STF at the Drosophila NMJ. However, we show clearly that a pV_{r}based facilitation mechanism (dual fusionsensor model) cannot account for STF in synapses with heterogeneous distances between release sites and Ca^{2+} channels. Interestingly, syt7 was also reported to function in vesicle priming and RRP replenishment (Liu et al., 2014; Schonn et al., 2008). Thus, future work will be necessary to investigate whether the function of syt7 in STF might take place by Ca^{2+}dependent inhibition of vesicle unpriming or release site activation.
Similar suggestions that facilitation results from a buildup of primed SVs during stimulus trains were made for the crayfish NMJ and mammalian synapses (Gustafsson et al., 2019; Pan and Zucker, 2009; Pulido and Marty, 2018). This is in line with our results, with facilitation arising from modulation of the number of primed SVs rather than pV_{r}. Our models are conceptually simple (e.g. all SVs are equally primed and distinguished only by distance to Ca^{2+ }channels, sometimes referred to as ‘positional priming’ Neher and Sakaba, 2008), and we improved conceptually on previous work by using estimated SV release site:Ca^{2+} channel distributions, stochastic simulations and comparison to variancemean relationships and we performed a systematic comparison of pV_{r} and primingbased models. It has not been clear whether increases in primed SVs are also required for pairedpulse facilitation, or only become relevant in the case of ‘tonic’ synapses that build up release during longer stimulus trains (frequency facilitation Neher and Brose, 2018). Pairedpulse facilitation is a more widespread phenomenon in synapses than frequency facilitation, and we show here for the case of Drosophila NMJ that it also seems to require primingbased mechanisms. Thus, Ca^{2+}dependent increases of the RRP during STP might be a general feature of chemical synapses.
Materials and methods
Fly husbandry, genotypes and handling
Request a detailed protocolFlies were kept under standard laboratory conditions as described previously (Sigrist et al., 2003) and reared on semidefined medium (Bloomington recipe) at 25°C, except for GCaMP6m and synapGCaMP6f flies which were kept at room temperature, and Ok6GAL4/+ (Figure 2, Figure 2—figure supplement 1, Figure 4 panel BE and G, Figure 6 panel DG and I, Figure 6—figure supplement 1, Figure 7D–G and I, Figure 7—figure supplement 1, Figure 7—figure supplement 3C–F,H) which were kept at 29°C (for detailed genotypes see below). For experiments both male and female 3^{rd} instar larvae were used. The following genotypes were used:
Figure 7 Ok6GAL4/+ (Ok6Gal4/II crossed to w[1118]; panel DG, (I). Figure 7—figure supplement 1: Ok6GAL4/+ (Ok6Gal4/II crossed to w[1118]). Figure 7—figure supplement 3: Ok6GAL4/+ (Ok6Gal4/II crossed to w[1118]; panel CF, (H).
The following stocks were used: Ok6GAL4/II (Aberle et al., 2002), UASUnc13AGFP/III (Böhme et al., 2016), elavGal4/I (Lin and Goodman, 1994). The following stock were obtained from the Bloomington Drosophila Stock Center: P{w[+mC]=MhcSynapGCaMP6f}3–5/III (Newman et al., 2017) and w[1118]; P{y[+t7.7] w[+mC]=20XUASIVSGCaMP6m}attP40. The following stock was obtained from Kyoto Stock Center: P84200/IV.
EM data acquisition and analysis
Request a detailed protocolSample preparation, EM image acquisition and the quantification of docked SV distances to the AZ center (center of the electron dense ‘Tbar’) are described in Böhme et al. (2016); ReddyAlla et al. (2017). The Rayleigh distributions were fit to the distances of docked SVs to the Tbar pedestal center, which had been collected in two EM datasets; analyses of these datasets were published in two previous studies, (ReddyAlla et al., 2017) for the histogram of distances depicted in Figure 1A and (Böhme et al., 2016) for the histogram of distances depicted in Figure 1—figure supplement 1A.
Derivation of the realistic docked SV distribution from EM measurements
Request a detailed protocolThe distances between Ca^{2+} channels and docked SVs in Drosophila NMJ obtained by EM was found to follow a Rayleigh distribution with best fit scale parameter σ = 76.51 nm (EM dataset 1) and σ = 74.07 nm (EM dataset 2). The fitting was performed with a MATLAB (MathWorks, version R2018b) function, raylfit, which uses maximum likelihood estimation. As these distances are found by EM of a crosssection of the active zone, we integrate this distribution around a circle to obtain the twodimensional distribution of SVs in the circular space around the active zone.
The Rayleigh distribution has the following probability density function (pdf):
The pdf of the SV distribution will then be a scaling of the following function
In order to find the pdf of the 2D SV distribution, we integrate $\widehat{g}$ to find the normalizing constant. By integration by parts we get
where the standard normal distribution was used in the last equality. Normalising (1) by this constant, we get the pdf of the distance distribution on a circular area in the active zone:
The SV distribution in simulations
Request a detailed protocolIn order to use the above SV distribution in simulations, we need to determine probabilities. g(x) is a generalized gamma distribution with $a=\sqrt{2}\cdot \sigma $, $p=2$, $d=3$. The generalized gamma distribution with a>0, p>0, d>0 has the following pdf:
and cumulative density function (cdf):
where $\gamma $ is the lower incomplete gamma function, and $\Gamma $ is the (regular) gamma function. Both of these functions are implemented in MATLAB (MathWorks, version R2018b), which easily allows us to draw numbers from them.
Thus, the SV distribution has the following cdf:
That is, given a uniformly distributed variable $q\in \left(\mathrm{0,1}\right)$, we can use inbuilt MATLAB functions to sample SV distances, d:
The implementation is as follows:
q = rand(1);
d = sqrt(2 * sigma ^ 2 * gammaincinv(q, 1.5));
Note that in MATLAB the inverse incomplete gamma function with parameter s is scaled by $\Gamma \left(s\right)$, which is why we input $q$ and not $q/\Gamma \left(1.5\right)$.
STED data acquisition and analysis
Request a detailed protocolSample preparation, Unc13A antibody staining, STED image acquisition and the isolation of single AZ images are described in Böhme et al. (2019) and in the following. Thirdinstar w[1118] larvae were put on a dissection plate with both ends fixed by fine pins. Larvae were then covered by 50 µl of icecold hemolymphlike saline solution (HL3, pH adjusted to 7.2 [Stewart et al., 1994]: 70 mM NaCl, 5 mM KCl, 20 mM MgCl_{2}, 10 mM NaHCO_{3}, 5 mM Trehalose, 115 mM DSaccharose, 5 mM HEPES). Using dissection scissors a small cut at the dorsal, posterior midline of the larva was made from where on the larvae was cut completely open along the dorsal midline until its anterior end. Subsequently, the epidermis was pinned down and slightly stretched and the internal organs and tissues removed. For the ‘STED dataset 2’ shown in Figure 1—figure supplement 1C,D, animals were then incubated in a HL3 solution containing 0.5% DMSO for 10 min (this served as a mock control for another experiment not shown in this paper using a pharmacological agent diluted in DMSO). The dissected samples were washed 3x with icecold HL3 and then fixed for 5 min with icecold methanol. After fixation, samples were briefly rinsed with HL3 and then blocked for 1 hr in 5% native goat serum (NGS; SigmaAldrich, MO, USA, S2007) diluted in phosphate buffered saline (Carl Roth Germany) with 0.05% TritonX100 (PBT). Subsequently dissected samples were incubated with primary antibodies (guineapig Unc13A 1:500; Böhme et al., 2016) diluted in 5% NGS in PBT overnight. Afterwards samples were washed 5x for 30 min with PBT and then incubated for 4 hr with fluorescencelabeled secondary antibodies (goat antiguinea pig STAR635 (1:100) diluted in 5% NGS in PBT. For secondary antibody production STAR635 fluorophore (Abberior, Germany) was coupled to respective IgGs (Dianova, Germany). Samples were then washed overnight in PBT and subsequently mounted in Mowiol (MaxPlanck Institute for Biophysical Chemistry, Group of Stefan Hell) on highprecision glass coverslips (Roth, Germany, LH24.1). Twocolor STED images were recorded on a custombuilt STEDmicroscope (Göttfert et al., 2017), which combined two pairs of excitation laser beams of 595 nm and 635 nm with one STED fiber laser beam at 775 nm. All STED images were acquired using Imspector Software (Max Planck Innovation GmbH, Germany). STED images were processed using a linear deconvolution function integrated into Imspector Software (Max Planck Innovation GmbH, Germany). Regularization parameter was 1e^{−11}. The point spread function (PSF for deconvolution was generated using a 2D Lorentz function with its halfwidth and halflength fitted to the halfwidth and halflength of each individual image. Single AZ images of ‘STED dataset 1’ (Figure 1E,F, Figure 1—figure supplement 1C,D) had previously been used for a different type of analysis defining AZ Unc13A cluster numbers; Wildtype in supplementary Figure 2a of Böhme et al. (2019). In this study here, we wanted to obtain the average Unc13A distribution from all AZs (no distinction of AZ types). To get an average image of the Unc13A AZ distribution, we used a set of hundreds of 51 × 51 pixel images with a pixel size of 10 × 10 nm. We identified Unc13A clusters in each image using the fluorescence peak detection procedure described in Böhme et al. (2019) using MATLAB (version 2016b). Peak detection was performed as follows: In each deconvolved 51 × 51 pixel image of an Unc13Astained AZ, a threshold of 25 gray values was applied below which no pixels were considered. Then, local maxima values were found by finding slope changes corresponding to peaks along pixel columns using the function diff. The same was done along rows for all column positions where peaks were found. The function intersect was then used to determine all pixel positions common in both columns and rows. A minimum distance of 50 nm between neighboring peaks was used to exclude the repeated detection of the same peak, and an edge of 10 nm around the image was excluded to prevent the detection of neighboring AZs. The center of mass of all peak x,ycoordinates found in a single image was then calculated as follows:
Here, n is the number of detected peaks, (P_{x}, P_{y}) represents the center of mass (x,y)coordinate, and x_{obs}(n) and y_{obs}(n) are the coordinates of the nth detected peak. The image was then shifted such that this position (P_{x},P_{y}) would fall into the center pixel of the 51 × 51 AZ image. For this, we calculated the required shift (d_{x} and d_{y}):
Here, imgsize(x,y) refers to the pixel dimensions of the image in both x and y dimensions. The required shift d_{x,y} was then applied to the image using imtranslate, which directly takes these shift values as an input. All shifted images were then averaged into a single compound average image of all AZs by taking the average of each individual pixel and linearly scaling the result in a range between 0 and 255. This resulted in a circular cloudy structure depicted in Figure 1E, Figure 1—figure supplement 1C. To obtain the distribution of fluorescence as a function of distance to the AZ center in the average picture, we determined the distance between the center of the image and the center of the pixel together with the fluorescence intensity in each pixel. The fluorescence intensity in each pixel was obtained by using the inbuilt MATLAB function ‘imread’, which outputs the intensities in a matrix with indexes corresponding to the pixel location in the picture. From the indexes (x_{p},y_{p}) of each pixel (of size 10 nm), the distance (in nm) to the center was calculated by the following formula:
We subtracted 26 from the pixel number, since the center pixel is the 26^{th} pixel in x and ydirection.These distances together with the intensity at each pixel provided the data for the histograms in Figure 1F and Figure 1—figure supplement 1D. The intensity values were normalized to the total amount of intensity making the yaxis of the histogram show percentage of the total amount of intensity.
Calculation of mean distance to four nearest neighbors (1–4NND)
Request a detailed protocolStage L3 larvae (n = 17; genotype: w[1118]; P{w[+mC]=MhcSynapGCaMP6f}3–5, Bloomington #67739) were fixed in icecold Methanol for 7 min and IHCstained for BRP (mouse antiNc82, 1:1000; secondary AB: goat antimouse Cy5 1:500). Confocal images of the preparations were taken and processed as described in ReddyAlla et al. (2017) for a different set of experiments not shown in this paper. Subsequently, the BRP channel was used to identify local fluorescence intensity maxima using the ImageJfunction ‘Find Maxima’ with a threshold setting between 10 and 20. The locations of maxima for each cell were then loaded into MATLAB (version 2016b) and the distances of each x,ycoordinate to all others were determined using the MATLAB function pdist2, resulting in a square matrix containing all possible interAZ distances. Each column of this matrix was then sorted in ascending order, and (as the distance of one AZ to itself is always 0) the mean of the 2^{nd} to 5^{th} smallest values across all AZs was determined and depicted as 1NND through 4NND in Figure 3A. The mean distance of the four nearest neighbouring AZs (1–4NND) was calculated in each AZ (gray circles in Figure 3A bottom right) and the mean across AZs was used for quantification of the simulation volume (see below).
Electrophysiological data acquisition and analysis
For both eEJC and mEJC (spontaneous release events,”miniature Excitatory Junctional Currents’) recordings, two electrode voltage clamp (TEVC) recordings were performed from muscle 6 NMJs of abdominal segments A2 and A3 as reported previously (Qin et al., 2005). Prior to recordings, the larvae were dissected in haemolymphlike solution without Ca^{2+} (HL3, pH adjusted to 7.2 Stewart et al., 1994: 70 mM NaCl, 5 mM KCl, 20 mM MgCl_{2}, 10 mM NaHCO_{3}, 5 mM Trehalose, 115 mM DSaccharose, 5 mM HEPES) on Sylgard (184, Dow Corning, Midland, MI, USA) and transferred into the recording chamber containing 2 ml of HL3 with CaCl_{2} (concentrations used in individual experiments described below). TEVC recordings were conducted at 21°C using sharp electrodes (borosilicate glass with filament, 0.86×1.5×80 nm, Science Products, Hofheim, Germany) with pipette resistances between 20–30 MΩ, which were pulled with a P97 micropipette puller (Sutter Instrument, CA, USA) and filled with 3 mM KCl. Signals were lowpass filtered at 5 KHz and sampled at 20 KHz. Data was obtained using a Digidata 1440A digitizer (Molecular devices, Sunnyvale, CA, USA), Clampex software (v10.6) and an Axoclamp 900A amplifier (Axon instruments, Union City, CA, USA) using Axoclamp software. Only cells with a resting membrane potential V_{m} below −50 mV, membrane resistances R_{m} above 4 MΩ and an absolute leak currents of less than 10 nA were included in the dataset.
eEJC recordings
Request a detailed protocoleEJC recordings were conducted at a membrane holding potential of −70 mV in TEVC mode. APs were evoked by giving 300 µs short depolarizing pulses (8 V) to respective innervating motoneuron axons using a suction electrode (pulled with DMZUniversal Puller (ZeitzInstruments GmbH, Germany) polished with the CPM2 microforge (ALA Scientific, NY, USA)) and a stimulator (S48, Grass Technologies, USA).
For experiments shown in Figure 2, individual cells were recorded at an initial extracellular CaCl_{2} concentration of 0.75 mM which was subsequently increased to 1.5 mM, 3 mM, 6 mM and 10 mM by exchanging and carefully mixing 1 ml of the bath solution with 1 ml HL3 of a higher CaCl_{2} concentration (total concentrations of exchange solutions: 2.25 mM, 4.5 mM, 9 mM, 14 mM), ultimately adding up to the desired CaCl_{2} concentration in the bath. At each titration step, cells were acclimated in the bath solution for 60 s and 10 repetitions of paired stimulating pulses (0.1 Hz, 10 ms interstimulus interval) were given. eEJC data shown in Figure 2—figure supplement 3 was obtained by recording Ok6Gal /+ and +/+ NMJs at 0.75 mM (Figure 2—figure supplement 3AD ) and 1.5 mM (Figure 2—figure supplement 3EH ) Ca^{2+}. A single test AP was given (followed by a 20 s intermission) and cells were stimulated once by two consecutive APs (10 ms interstimulus interval). In Figure 2—figure supplement 3B, D, E, and G, eEJC_{1} and PPR averages are shown ± the estimated singlecell SD .
eEJC data was analyzed with our own custombuilt MATLAB script (provided with the source data file, Figure 2—source data 1). After stimulation artifact removal, the eEJC_{1} amplitude was determined as the minimum current value within 10 ms from the time of stimulation. To account for the decay only being partial before the second stimulus, we fitted a single exponential function to the eEJC decay from the time point of 90% of the amplitude to the time point of the second stimulus. The eEJC_{2} amplitude was determined as the difference between the minimum after the second stimulus and the value of the fitted exponential at the time point of the second minimum (see insert in Figure 2C and Figure 2—figure supplement 1A). For analysis shown in Figure 2, the first stimulation per Ca^{2+} concentration was excluded, as we noticed that the first trial often gave first eEJC responses that were higher than in the following trials. This may reflect the presence of a slow reaction by which SVs can be primed with an even higher release probability (possibly due to the ‘superpriming’ described at the murine Calyx of Held synapse Lee et al., 2013). However, as the var/mean analysis requires the existence of an equilibrium inbetween stimuli which appears to have been reached between all of the succeeding stimuli, we decided to use only those for our analysis. For eEJC_{1} amplitudes the average over all measurements and all cells (6 cells, nine measurements each) was calculated (Figure 2B). The PPR was calculated by dividing the second amplitude by the first throughout trials and averaging over all measurements and all cells (Figure 2D). In each cell, the variance of eEJC_{1} and PPR was estimated (nine stimulations per Ca^{2+} concentration) and the average variance (averaged across cells) was calculated at each extracellular Ca^{2+} concentration. The error bars in Figure 2B,D are the SD (across all animals) at each extracellular Ca^{2+} concentration. In Figure 2F the eEJC_{1} averages and variances are ± SEM. A parabola with intersect y = 0 was fitted using the function polyfitZero (version 1.3.0.0 from MathWorks file exchange) in MATLAB. (Var = q*II^{2}/N, q being the quantal size, I the mean eEJC_{1} amplitude and N number of release sites) (Clements and Silver, 2000).
mEJC recordings
Request a detailed protocolmEJC data was obtained from a separate set of experiments where mEJCs were recorded for 60 s in TEVC mode at 1.5 mM extracellular Ca^{2+} and a holding potential of −80 mV for easier identification of miniature events. Because different holding potentials were used (−80 mV here compared to −70 mV for the data shown in Figure 2) it must be pointed out that these recordings were only used to determine the shape of the response for later convolution with SV fusion events predicted by the model (the mEJC amplitude was adjusted based on the variancemean data collected at 70 mV, see below). For this, the average mEJC traces from five different cells were aligned to 50% of the rise and averaged. We then fitted the following formula to the data:
${t}_{0}$ is the onset, $A$ is the full amplitude (if there was no decay), $B$ is the fraction of the fast decay, and ${\tau}_{r},{\tau}_{df},{\tau}_{ds}$ are the time constants of the rise, fast decay, and slow decay respectively.
The best fit was
and is plotted together with the average experimental mini trace in Figure 2—figure supplement 1B. Note that ${t}_{0}$ is a time delay when this mEJC is implemented in the simulation and is therefore arbitrary. B is very small making the decay close to a single exponential. The maximum of this function is ~0.7 nA. However, as mentioned above, this function was rescaled to a value of 0.6 nA to match the mEJC amplitudes of the experiments conducted with a holding potential of 70 mV, that is the size of a single quantal event, q=0.6 nA, estimated from the variancemean analysis (see Figure 2F).
Presynaptic GCaMP recordings and analysis
Request a detailed protocolBecause the presynaptic terminals of the Drosophila larval NMJ are not readily accessible to electrical recordings of Ca^{2+} currents, the saturation behaviour of Ca^{2+} influx as a function of extracellular Ca^{2+} concentrations was measured. We did so by engaging the fluorescent Ca^{2+} indicator GCaMP6m (Genotype: w[1118]; P{y[+t7.7] w[+mC]=20XUASIVSGCaMP6m}attP40, Flybase ID: FBti0151346), which we expressed presynaptically using OK6Gal4 as a motoneuronspecific driver. Third instar larvae heterozygously expressing the indicator were used in experiments as follows. Dissection took place in Ca^{2+}free, standard hemolymphlike solution HL3 (in mM: NaCl 70, KCl 5, MgCl_{2} 20, NaHCO_{3} 10, Trehalose 5, Sucrose 115, HEPES 5, pH adjusted to 7.2) (Stewart et al., 1994). After dissection on a Sylgard184 (DowCorning) block, larvae were transferred to the recording chamber containing HL3 at varying CaCl_{2} concentrations (see below). The efferent motoneuron axons were sucked into a polished glass electrode containing a chlorided silverwire, which could be controlled via a mechanical micromanipulator (Narishige NMN25) and was connected to a pipette holder (PPH1PBNC, NPI electronics) via a patch electrode holder (NPI electronics), and connected to an S48 stimulator (Grass Technologies). Larvae were then recorded using a whitelight source (Sutter DG4, Sutter Instruments) and a GFP filter set with a Hamamatsu OrcaFlash 4.0v2 sCMOS (Hahamatsu Photonics) with a framerate of 20 Hz (50 ms exposure) controlled by µManager software (version 1.4.20, https://micromanager.org) on an upright microscope (Olympus BX51WI) with a 60x waterimmersion objective (Olympus LUMFL 60 × 1.10 w). Muscle 4 1b NMJs in abdominal segments 2 to 4 were used for imaging. Imaging was conducted over 10 s, and at 5 s, 20 stimuli were applied to the nerve at 20 Hz in 300µs 7V depolarization steps. This procedure was begun in the lowest Ca^{2+} concentration (0.75 mM) and then repeated in the same larva at increasing Ca^{2+} concentrations (in mM 1.5, 3, 6) by exchanging the extracellular solution. To achieve a situation with no Ca^{2+} influx, a final recording was conducted where the bath contained HL3 without CaCl_{2} and instead 8.3 mM EGTA (this solution was made by diluting 2.5 ml of a 50 mM stock solution in H_{2}O in 12.5 ml of HL3, resulting in a pH of 8.0). Because this results in a slight dilution (16%) of the components in the HL3, the same dilution was performed for the above described Ca^{2+}containing solutions by adding 2.5 ml H_{2}O to 12.5 ml of HL3 before CaCl_{2} was added at above mentioned concentrations.
Analysis of 5 Drosophila 3^{rd} instar Larvae was done after automated stabilization of x,ymovement in the recordings (8bit multipage .TIFstacks, converted from 16 bit) as described previously (ReddyAlla et al., 2017), manually selecting a ROI around the basal fluorescent GCaMP signal, and reading out the integrated density (the sum of all pixel grey values) of the whole region over time. Background fluorescence was measured in a region of the same size and shape outside of the NMJ and subtracted (framewise) from the signal, separately for each single recording. The quantification was then performed individually for each Ca^{2+} concentration, by subtracting the fluorescence 250 ms before the stimulation (F_{t=4.75s}) from the maximum fluorescence of the trace (F_{max}), yielding the change in fluorescence dF:
This was repeated for each cell and a Hill fit was performed on the individual values using Prism (version 6.07, GraphPad Software Inc):
In the above equation, F_{end} is the asymptotic plateau of the fluorescence increase. Furthermore, [Ca^{2+}]_{ext} is the extracellular Ca^{2+} concentration. K_{M,fluo} (best fit value: 2.679 mM) is the concentration of extracellular Ca^{2+} at which fluorescence was half of F_{end}. The exponent m indicates a cooperative effect of the extracellular Ca^{2+} concentration on the fluorescence increase, which was constrained to a value of 2.43 (unitless) based on the described Ca^{2+ }cooperativity of GCaMP6m (Barnett et al., 2017). However, constraining this value only had a modest effect on the estimate of K_{M,fluo} as leaving it as a free parameter yielded similar values for K_{M,fluo} (3.054 mM) and m (1.887). The constant C added at the end of Equation 3 allowed the baseline fluorescence to be different from zero. Results and best fit are summarized in (Figure 3—figure supplement 1).
Proof that stochastic simulation of release is needed for PPR estimation
Request a detailed protocolWe here prove that stochastic simulations of neurotransmitter release provide a different average PPR value than the PPR value estimated in deterministic simulations. In the following, the stochastic variables A_{1} and A_{2} represent the amplitudes of the first and second release, respectively, capital ‘E’ denotes the mean of a stochastic variable (e.g. EA_{1}), and a_{1} and a_{2} represent the amplitudes of the first and second release in the deterministic simulations. In all cases of parameter sets that we tried, the average amplitudes from the stochastic simulations with 1000 repetitions differed < 0.5 nA from the deterministically determined amplitudes. Thus, we can assume that EA_{1 }= a_{1} and likewise for the second release.
In deterministic simulations, the estimate of the PPR is
On the other hand, stochastic simulations yield a sample of different PPR values, since repetitions of the simulation routine yield release varying from trial to trial. In that case, the estimated PPR is
This resembles the way the PPR is estimated in experiments.
Using Jensen’s Inequality and the fact that the function f(x)=1/x is strictly convex, we get
Applying this to (4) we get
Thus, the average stochastically simulated PPR do not necessarily converge to the deterministic estimate with increasing repetitions (note that in general it is true that the mean of a nonlinear function of two random variables is not equal to the nonlinear function evaluated in the means). An example is shown in Figure 4—figure supplement 1, where the singlesensor model was simulated with varying amounts of Ca^{2+} influx (by varying Q_{max}). The most left blue point, for example, is significantly higher than the deterministic estimate (p=4e16, onesample ttest). This motivates the use of stochastic simulations for correct estimation of the PPR.
Simulation flow
Request a detailed protocolAll MATLAB procedures for simulation of the models can be found in Source code 1.
All simulations (deterministic and stochastic, see below) consisted of the same four basic steps, which we describe in detail here.
Given a set of parameters, we first ran deterministic Ca^{2+} simulations in space and time in the presynapse at the desired extracellular Ca^{2+ }concentrations.
A set of SV distances was drawn from the generalized gamma distribution. The set of SV distances provided the points at which to read the intracellular Ca^{2+} concentrations for the exocytosis simulation.
The simulation of the models for Ca^{2+} binding and exocytosis was performed for each SV position with the Ca^{2+} transients giving rise to the changing reaction rates.
The outcome of the exocytosis simulation were convolved with a mEJC which yielded the eEJC.
For each new set of parameters, steps 1–4 were repeated. For stochastic simulations, steps 2–4 were repeated 1000 times except for the parameter exploration in Figures 6J–K and 7J–K, where we ran 200 repetitions per parameter set. The many repetitions allowed a good estimate of both mean and variance of the models. In all cases, the mean amplitudes from the stochastic simulations with 1000 repetitions differed < 0.5 nA from the deterministically determined amplitudes.
Ca^{2+} simulation
Request a detailed protocolSimulation of Ca^{2+} signals in the presynapse was performed with the program CalC version 6.8.6 developed and maintained by Victor Matveev (Matveev et al., 2002). After this work was initiated, a bug affecting simulations of multiple Ca^{2+} channels in the same topology was found and a new version of CalC was released. This update had no effect on the simulations used in this study.
Intracellular Ca^{2+} concentrations were simulated in space and time in a cylindershaped volume. The cylinder allowed us to assume spatial symmetry which reduced simulation time significantly. Borders of the simulation volume were assumed to be reflective to mimic diffusion of Ca^{2+} from adjacent AZs (Meinrenken et al., 2002) and a volumedistributed uptake mechanism was assumed.
From measurements of the distance between an AZ and its four nearest neighbors (Figure 3A) we estimated the distance between centers of active zones to be 1.106 µm, leading to the assumption that the AZ spans a square on the membrane with area of 1.223 µm^{2}. In order for the cylindrical simulation volume to cover an area of the same size, the radius was set to 0.624 µm. The height of the simulation volume was set to 1 µm making the simulation volume 1.223 µm^{3}. Increasing the height further had no effect on the Ca^{2+} transients.
The total amount of charge flowing into the cell was assumed to relate to extracellular Ca^{2+} in a MichaelisMentenlike way (as previously described by Schneggenburger et al., 1999; Trommershäuser et al., 2003) such that
K_{M,current} was set to the value of 2.679 mM as determined for K_{M,fluo} in the GCaMP6m experiments (see above). Q_{max} was fitted during the optimizations of the models.
We simulated a 10 ms paired pulse stimulus initiated after 0.5 ms of simulation. The Ca^{2+} currents for the two stimuli were simulated for 3 ms each and assumed to be Gaussian with FWHM = 360 µs and peak 1.5 ms after initiation. That is:
with $\sigma =\text{}\frac{0.360}{2\sqrt{2\cdot \mathrm{l}\mathrm{n}\left(2\right)}}=\text{}\mathrm{0.153.}$
The CalC simulation output were data files that contained the spatiotemporal intracellular Ca^{2+} profile at the height of 10 nm from the plasma membrane. In exocytosis simulations, these concentrations were interpolated at the SV distances in the x,yplane and at time points with MATLAB’s builtin interpolate functions when computing the reaction rates of the system at a given time point.
The resting Ca^{2+} concentration was assumed to relate to the extracellular Ca^{2+} concentration in a similar way as during stimulation, such that
with ${\left[C{a}^{2+}\right]}_{max}=190nM$
For designation and value of Ca^{2+} parameters, see Table 1.
SV distribution drawing
Request a detailed protocolIn all simulations we had to determine where to place release site. This was done by using the cdf of the SV distance distribution derived above (Equation 2).
For deterministic simulations, which were used in the fitting routine of the models (see below), the unit interval was divided into 180 bins of the form
The midpoints were the percentiles giving rise to distances at which we read the Ca^{2+} simulation. This approach provided an approximation of the SV distribution. In accordance with our assumption that the AZs work in parallel the 180 distances gave rise to 180 independent different systems of ODEs with 1/180 of the total amount of SVs in each system. The results were then added together as a good approximation of the mean of the stochastic simulations with random SV distance drawings.
In each run of the stochastic simulations, we drew n random numbers from the unit interval, n being the number of SVs, and computed the distances based on the formula derived above.
Rate equations of the simulated models
Request a detailed protocolThe models are summarized in Figures 4A, 6A and 7A, and Figure 7—figure supplement 3A,B. In the following equations the singlesensor, dual fusionsensor, and unpriming models are all described. The site activation model is a combination of the equations for the singlesensor model and the site activation equations described below. The denotes terms that are unique to the dual fusionsensor model, indicates unpriming, which is unique to the unpriming model. Parameters are described below. For designation and value of parameters, see Tables 2,3.
Rate equations of the singlesensor model, dual fusionsensor model and unpriming model:
In the singlesensor and site activation models, k_{2} = k_{2}=u = 0, and s = 1. This excludes all reactions exclusive for the dual fusionsensor and unpriming models. Similarly, u = 0 in the dual fusionsensor model and k_{2} = k_{2}=0 and s = 1 in the unpriming model.
[R(n,m)] denotes the Ca^{2+} binding state of a SV with n Ca^{2+} ions bound to the first sensor and m Ca^{2+} ions bound to the second fusion sensor. Note that in the singlesensor, site activation and unpriming models, m is always zero (since there is no second fusion sensor), and the states are denoted with a single number in Figures 4A and 6A and Figure 7—figure supplement 3. [F] counts the cumulative number of fused SVs. [P0] is not shown in the figures, but are part of the equations denoting the number of empty sites. That is, in the singlesensor and unpriming models $r=1\frac{{\left[C{a}^{2+}\right]}^{n}}{{\left[C{a}^{2+}\right]}^{n}+{K}_{M,unprim}^{n}}$ has a positive part equal to $\frac{d\left[P0\right]}{dt}$ and a negative part equal to the rate of replenishment. In the dual fusionsensor model, there are three states of empty sites, [P0], [P1], [P2]. These corresponded to the different states of Ca^{2+} binding to the second fusion sensor of the empty sites since we assumed the second sensor to be located on the plasma membrane. Note that these equations describe the second sensor with cooperativity 2, which is described in Results. We also optimized cooperativities 3, 4, and 5. The equations can easily be extended to these cases, since the rate equations of the second fusion sensor are of the same form as for the first sensor. In the unpriming model (Figure 7A) we assumed unpriming to take place from state [R(0)] with a Ca^{2+}dependent rate.
For the individual reactions, we can express the rates of Ca^{2+} (un)binding, fusion, and replenishment of a single SV in a more general form. This is useful in the stochastic simulation method introduced later. In the following, we denote the general form of the rate for each possible reaction in the models described above.
The expressions in brackets denote the states involved in the reaction.
with n_{max} and m_{max} denoting the cooperativity of the first and second fusion sensors, respectively. Equations in line 3 and 4 in (7) were only nonzero in the dual fusionsensor model.
Rate equation of the site activation model
Request a detailed protocolIn the site activation model (Figure 7—figure supplement 3), all reactions regarding Ca^{2+} (un)binding and replenishment was as in the onesensor model. In addition we assumed a mechanism acting on the release sites independently of the Ca^{2+} binding of the SV. All sites regardless of the SV status were either activated (A state) or not (D or I states). This mechanism is proposed as a facilitation mechanism, which necessitates its primary effect to be on the second stimulus rather than the first. We were therefore forced to implement the D state, which is a temporary ‘delay’ state making sure the mechanism does not increase first release. The changing of [A] and [I] states at 0.75 and 10 mM extracellular Ca^{2+} are shown in (Figure 7—figure supplement 3I).
The site activation mechanism has the following rate equations:
where $\alpha ,\text{}\beta ,\text{}\delta ,\text{}\gamma 0$ are rate parameters.
The deterministic implementation of the site activation model included 3 sets of ODEs, one for each state in the site activation model. Each set consisted of the equations of the onesensor model as well as transitions between states of equal Ca^{2+} binding in the 3 sets of ODEs (e.g. from R(0,D) to R(0,A)) (Figure 7—figure supplement 3B).
In the stochastic simulations the site activation rates were included in the propensity vector like any other reaction. Whenever a site activation reaction occurred, a release site vector consisting of n_{sites} elements was updated. For each site, the fusion rate was multiplied by 0, when the site state was I or D.
Steadystate estimation
Request a detailed protocolPrior to simulation, the Ca^{2+} binding states of all SVs were assumed to be in equilibrium. We can determine the steady state iteratively by setting
This can be reduced to the noniterative expression:
Note that for n = 0, the first parenthesis is 1, while m = 0 implies that the second parenthesis is 1, making this solution valid also in the absence of a second fusionsensor. We ignored the very small fusion rate. In the steadystate of the unpriming model, the number of SVs in [R(0,0)] must furthermore be in equilibrium with the number of empty states:
After finding this steadystate, the solution is scaled to match the desired number of SVs by multiplying all states with a constant, such that the sum of all [R(n,m)] and [P] equals the number of SVs. The steadystate of the site activation was determined before simulation by calculating the fraction of states being in [A], [D], or [I]. This was done by calculating
and normalizing to sum to 1. This determined the steady state fraction of activation of sites. In the stochastic simulations, the SVs were randomly assigned initial states according to the probabilities of the different states in the steadystate.
Deterministic exocytosis simulation
Request a detailed protocolAll deterministic exocytosis simulations of the above equations were carried out with the inbuilt MATLAB ODE solver ode15s.
Stochastic exocytosis simulation
Request a detailed protocolAll stochastic exocytosis simulations as well as simulation data handling were carried out in MATLAB with customwritten scripts (included in Source code 1). For the simulation itself we used a modified version of the Gillespie Algorithm (Gillespie, 2007), which included a minimal time step since reaction rates change quickly with the changing intracellular Ca^{2+} concentration. The minimal step was µ = 1e6 s. In the algorithm, the time from the current simulation time point, t, until the next reaction,τ, is determined, the reaction is carried out and the new simulation time point is set to t+τ. Whenever the simulation yielded τ>µ, the simulation time point was set to t+µ, no reaction was carried out and the propensities of the model were updated at the new time point. This is a valid method of obtaining a better estimate because the waiting time until next reaction is exponentially distributed.
The implementation of the algorithm takes advantage of the general form of the rate equations in (7). Instead of calculating matrices of states and reaction rates, we have a vector, V, of length n_{sites}, where each element represents the status of one SV/site. The SV state of a docked SV on the k^{th} site in state [R(n,m)] is denoted by the twodigit number
If the site was empty (due to initial submaximal priming or SV fusion) we assigned ${V}_{k}=100$.
Using Equation 7, the rates of any primed SV are
The sum of these rates of all SVs yield the summed propensities of the system, a_{0}, which is the basis of the calculation of τ, whereas the cumulative sum is used for determination of which SV undergoes a reaction (Gillespie, 2007). When a SV undergoes a reaction, we find the index of the reaction occurring, j, by using the cumulative sum of r_{k} in the same way as in the standard implementation of the Gillespie Algorithm (Gillespie, 2007). Putting $\hat{j}=j3$ allows us to easily update the status of the SV, since
In parallel with this a vector of fusions is updated, such that at every time point, the next element in the fusion vector is set to 1 if a fusion took place, and 0 else.
Parallel computing
Request a detailed protocolMany repetitions of time consuming stochastic simulations had to be performed, and many sets of ODEs were solved for each choice of parameters. Therefore, simulations were carried out on the computer grid on The Bioinformatics Center, University of Copenhagen. This allowed running repetitions in parallel with MATLAB’s Parallel Computing toolbox using between 5 and 100 cores depending on the simulation job.
Calculating the postsynaptic response
Request a detailed protocolIn order to calculate the eEJC, we needed a vector of the SV fusions at different time points. Both deterministic and stochastic simulations yielded the vectors time_outcome and fuse_outcome, which is a pair of vectors of the same length but with changing time steps. For the sampling we generated a time vector, time_sample, with a fixed time step of 1 µs. From here, the determining of the SV fusion times differ between deterministic and stochastic simulations.
In the deterministic simulations, we simulated a sample of distances, bins, as described earlier. Each bin gave rise to a set of ODEs, which could be simulated independently, and the fuse_outcome is continuously changing based on the rates. In MATLAB the interpolation for bin k was done as follows:
fuse_interp_{k} contained the cumulative fused SVs over time in a single bin sampled at the time points of the vector time_sample. These were summed to find the total number of fused SVs:
Therefore the SVs fused per time step were be the difference between neighboring values in the fuse_interp vector:
This vector was the basis for the computation of the eEJC.
In the stochastic simulations, the fuse_outcome vector contains discrete SV fusions at certain time points. We therefore sample the SV fusions by assigning them to the nearest time points on the time_sample vector. That is, each fusion time was rounded to the nearest microsecond, thereby giving rise to the fusion_vec, which in the stochastic case contained whole numbers of SV fusions at different time points.
In both deterministic and stochastic simulations the mEJC was generated as a vector, mEJC_vec, with the same time step as the time_sample and fuse_vec. This allows us to calculate the eEJC with MATLAB’s convolve function, conv, such that
where fusion_vec is a vector with the same time step, each element being the number of SV fusions at each time point.
Analysis of simulated eEJCs
Request a detailed protocolThe eEJC_{1} amplitude was determined as the minimum current of the eEJC within the time interval (0,10) ms. Similar to the analysis of experimental eEJC data, we fitted an exponential function to the decay for estimation of the base value for the second response (see Figure 2—figure supplement 1A). The eEJC_{2} amplitude was the difference between the second local minimum and the fitted exponential function extrapolated to the time point of the second local minimum (as described for the analysis of electrophysiology experiments).
Fitting routine
Request a detailed protocolBecause deterministic simulations cannot predict PPR values (due to Jensen’s inequality, see above), but stochastic simulations cannot be fitted to data, we first ran deterministic simulations comparing the simulated first and second absolute eEJC amplitudes to the experimental amplitudes (not the PPR, see Materials and methods). Afterwards we ran stochastic simulations with the optimised parameters in order to compare PPRs and variances to experimental results. To determine the optimal parameters for the deterministic simulations at the five experimental extracellular Ca^{2+} concentrations, the models were fitted to the two peak amplitudes, eEJC_{1} and eEJC_{2}, by minimizing the following cost value:
where we sum over the five different experimental Ca^{2+} concentrations. Note that in deterministic simulations, eEJC_{1} and eEJC_{2} amplitudes are precise estimates of average amplitudes in stochastic simulations allowing us to do deterministic optimizations.
When fitting the models, we used the inbuilt MATLAB function fminsearch, which uses the NelderMead Simplex Search, to minimize the above cost function. The cost calculation in each iteration was a twostep process taking advantage of the fact that the total number of SVs scales the eEJC_{1} and eEJC_{2} values in the deterministic simulations. For each choice of parameters the simulation was run with 180 sites (the initial number of sites is arbitrary, but matched the number of distance bins), and the optimal number of sites were determined afterwards. Thus, a given set of parameters gave rise to amplitudes eEJC_{1,init} and eEJC_{2,init} from simulations with 180 sites. After that we determined $eEJC=conv(fuse\_vec,mEPSC\_vec)$ such that $cost(eEJ{C}_{1,sim},eEJ{C}_{2,sim})=\sum _{k=1}^{5}\left(\frac{{\left(eEJ{C}_{1,sim,k}eEJ{C}_{1,exp,k}\right)}^{2}}{eEJ{C}_{1,exp,k}}+\frac{{\left(eEJ{C}_{2,sim,k}eEJ{C}_{2,exp,k}\right)}^{2}}{eEJ{C}_{2,exp,k}}\right)$ was minimized. The number of sites in the given iteration was therefore 180⋅c_{sites} and the cost of that particular iteration was
In this way the optimization algorithm did not have to include n_{sites} in the parameter search algorithm, which reduced the number of iterations significantly.
In the stochastic simulations, the number of SVs was set to 180⋅c_{sites} rounded to nearest integer.
References
 1
 2
 3
 4
 5
 6
 7
 8

9
Rapid active zone remodeling consolidates presynaptic potentiationNature Communications 10:1085.https://doi.org/10.1038/s41467019089776
 10

11
Control of synaptic strength and timing by the releasesite Ca2+ signalNature Neuroscience 8:426–434.https://doi.org/10.1038/nn1417

12
Facilitation of presynaptic calcium currents in the rat brainstemThe Journal of Physiology 513 ( Pt 1:149–155.https://doi.org/10.1111/j.14697793.1998.149by.x
 13

14
Fife organizes synaptic vesicles and calcium channels for highprobability neurotransmitter releaseThe Journal of Cell Biology 216:231–246.https://doi.org/10.1083/jcb.201601098
 15
 16
 17

18
Single Ltype calcium channel conductance with physiological levels of calcium in chick ciliary ganglion neuronsThe Journal of Physiology 496 ( Pt 1:59–68.https://doi.org/10.1113/jphysiol.1996.sp021665

19
Unveiling synaptic plasticity: a new graphical and analytical approachTrends in Neurosciences 23:105–113.https://doi.org/10.1016/S01662236(99)015209
 20

21
Cooperative action a calcium ions in transmitter release at the neuromuscular junctionThe Journal of Physiology 193:419–432.https://doi.org/10.1113/jphysiol.1967.sp008367

22
Nanodomain coupling between Ca2+ channels and sensors of exocytosis at fast mammalian synapsesNature Reviews Neuroscience 13:7–21.https://doi.org/10.1038/nrn3125
 23

24
Shortterm forms of presynaptic plasticityCurrent Opinion in Neurobiology 21:269–274.https://doi.org/10.1016/j.conb.2011.02.003
 25

26
Maturation of active zone assembly by Drosophila bruchpilotThe Journal of Cell Biology 186:129–145.https://doi.org/10.1083/jcb.200812150

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

30
The first 100 nm inside the Presynaptic terminal where calcium diffusion triggers vesicular releaseFrontiers in Synaptic Neuroscience 10:23.https://doi.org/10.3389/fnsyn.2018.00023
 31
 32
 33
 34

35
A mathematical model of synaptotagmin 7 revealing functional importance of shortterm synaptic plasticityNeural Regeneration Research 14:621–631.https://doi.org/10.4103/16735374.247466
 36
 37
 38
 39
 40

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

42
The role of calcium in neuromuscular facilitationThe Journal of Physiology 195:481–492.https://doi.org/10.1113/jphysiol.1968.sp008469
 43
 44
 45
 46
 47
 48

49
Synaptotagmin controls and modulates synapticvesicle fusion in a ca(2+)dependent mannerTrends in Neurosciences 18:177–183.https://doi.org/10.1016/01662236(95)938988
 50
 51

52
New insights into shortterm synaptic facilitation at the frog neuromuscular junctionJournal of Neurophysiology 113:71–87.https://doi.org/10.1152/jn.00198.2014
 53
 54

55
The bruchpilot cytomatrix determines the size of the readily releasable pool of synaptic vesiclesThe Journal of Cell Biology 202:667–683.https://doi.org/10.1083/jcb.201301072

56
New and corrected simulations of synaptic facilitationBiophysical Journal 83:1368–1373.https://doi.org/10.1016/S00063495(02)739076
 57

58
Residual bound Ca2+ can account for the effects of Ca2+ buffers on synaptic facilitationJournal of Neurophysiology 96:3389–3397.https://doi.org/10.1152/jn.00101.2006

59
Calcium secretion coupling at Calyx of held governed by nonuniform channelvesicle topographyThe Journal of Neuroscience 22:1648–1667.https://doi.org/10.1523/JNEUROSCI.220501648.2002
 60
 61

62
RIM controls homeostatic plasticity through modulation of the readilyreleasable vesicle poolJournal of Neuroscience 32:16574–16585.https://doi.org/10.1523/JNEUROSCI.098112.2012
 63
 64
 65
 66
 67
 68
 69

70
GAL4 drivers specific for type ib and type is motor neurons inDrosophilaG3: Genes, Genomes, Genetics 9:453–462.https://doi.org/10.1534/g3.118.200809
 71
 72

73
Quantal fluctuations in central mammalian synapses: functional role of vesicular docking sitesPhysiological Reviews 97:1403–1430.https://doi.org/10.1152/physrev.00032.2016

74
A twostep docking site model predicting different shortterm synaptic plasticity patternsThe Journal of General Physiology 150:1107–1124.https://doi.org/10.1085/jgp.201812072
 75
 76

77
The membrane fusion enigma: snares, Sec1/Munc18 proteins, and their accomplicesguilty as charged?Annual Review of Cell and Developmental Biology 28:279–308.https://doi.org/10.1146/annurevcellbio101011155818

78
Synaptic weight set by Munc131 supramolecular assembliesNature Neuroscience 21:41–49.https://doi.org/10.1038/s4159301700419
 79
 80
 81
 82
 83

84
Ca2+ from one or two channels controls fusion of a single vesicle at the frog neuromuscular junctionJournal of Neuroscience 26:13240–13249.https://doi.org/10.1523/JNEUROSCI.141806.2006

85
Experiencedependent strengthening of Drosophila neuromuscular junctionsThe Journal of Neuroscience 23:6546–6556.https://doi.org/10.1523/JNEUROSCI.231606546.2003
 86

87
The Nanophysiology of Fast Transmitter ReleaseTrends in Neurosciences 39:183–197.https://doi.org/10.1016/j.tins.2016.01.005
 88

89
Improved stability of Drosophila larval neuromuscular preparations in haemolymphlike physiological solutionsJournal of Comparative Physiology A 175:179–191.https://doi.org/10.1007/BF00215114
 90

91
A molecular machine for neurotransmitter release: synaptotagmin and beyondNature Medicine 19:1227–1231.https://doi.org/10.1038/nm.3338
 92
 93
 94

95
Effects of mobile buffers on facilitation: experimental and computational studiesBiophysical Journal 78:2735–2751.https://doi.org/10.1016/s00063495(00)768196
 96
 97
 98

99
Simple stochastic models for the release of quanta of transmitter from a nerve terminalAustralian Journal of Statistics 8:53–63.https://doi.org/10.1111/j.1467842X.1966.tb00164.x
 100
 101

102
Shortterm plasticity at the Calyx of heldNature Reviews Neuroscience 3:53–64.https://doi.org/10.1038/nrn705
 103
 104

105
Multiple Ca2+ sensors in secretion: teammates, competitors or autocrats?Trends in Neurosciences 34:487–497.https://doi.org/10.1016/j.tins.2011.07.003
 106
 107
 108
 109
 110

111
Rapid active zone remodeling during synaptic plasticityJournal of Neuroscience 31:6041–6052.https://doi.org/10.1523/JNEUROSCI.669810.2011
 112
 113
 114
 115

116
Is synaptotagmin the calcium sensor?Current Opinion in Neurobiology 13:315–323.https://doi.org/10.1016/S09594388(03)000631

117
Changes in the statistics of transmitter release during facilitationThe Journal of Physiology 229:787–810.https://doi.org/10.1113/jphysiol.1973.sp010167

118
Shortterm synaptic plasticityAnnual Review of Physiology 64:355–405.https://doi.org/10.1146/annurev.physiol.64.092501.114547
Decision letter

Mark CW van RossumReviewing Editor; University of Nottingham, United Kingdom

Ronald L CalabreseSenior Editor; Emory University, United States

Victor MatveevReviewer; NJIT
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Acceptance summary:
The paper presents a thorough combination of EM and computational modelling to explain the short term plasticity of synaptic release.
Decision letter after peer review:
Thank you for submitting your article "Rapid regulation of vesicle priming explains synaptic facilitation despite variable vesicle:Ca^{2+} channel distances" for consideration by eLife. Your article has been reviewed by three peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by Ronald Calabrese as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Victor Matveev (Reviewer #1).
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission. While eLife aims to prevent long review cycles and substantial additional work, the reviewers have identified quite a number of distinct issues. In subsequent discussions these were all deemed relevant to support the study's conclusion. However, a number of issues can be addressed by providing a deeper discussion of the limitations and assumptions of the study.
Summary:
This work combines experimental recordings with mathematical modeling to investigate mechanisms of Ca^{2+}exocytosis coupling and facilitation (STF) at Drosophila NMJ. The distinguishing feature of this study is that vesiclechannel distances were carefully mapped, allowing the authors to investigate the implications of nonuniform channelvesicle distance on NT release at this NMJ. The main conclusion of this work is that many traditional models cannot account for the STF observed at this NMJ, given the observed distribution of channelvesicle distances. The reason is that traditional models explain STF through changes in the vesicle release probability, but the latter is very high close to the Ca^{2+} channel, and very low far from the Ca^{2+} channel. Therefore, very few vesicles would feel the effect of small Ca^{2+} accumulation. The authors propose and test an alternative model which explains the observed properties of neurotransmitter release at this terminal. This model posits that vesicle undergo priming prior to becoming fusionready. Although similar STF models have been recently considered by others (as cited in this work), the presented model is more comprehensive, implements important stochastic effects, and proposes inhibition of depriming rather than acceleration of priming as a key target of Ca^{2+} action. An alternative "site activation" model is also presented, and there is some parameter sensitivity analysis included.
The hypothesis presented for STF and Caexocytosis coupling is in general interesting, and the model comparisons presented contribute to the understanding of this problem.
Essential revisions:
1.1) Equation 2: The whole Ca^{2+} current calibration is very hard to follow, and is my main concern regarding methods description. In Equation 2, shouldn't K_{M} be raised to power "h"? Also, this expression gives the relationship between fluorescence and intracellular Ca^{2+}, or was the fluorescence recorded with extracellular GCaMP? Further, why is K_{M} determined from GCaMP6 calibration? Isn't K_{M} an innate property of the Ca channel conductivity (and possibly Cadependent inactivation), having nothing to do with the indicator dye? The entire Ca^{2+} current calibration has to be explained much more clearly, I could not follow it.
1.2) Unless I misunderstood the modeling/methods, the described stochastic simulation process, while technically correct, appears needlessly complicated and computationally expensive. Any Markov Chain with a single absorbing state and deterministic propensities ([Ca] is deterministic here) is described exactly by its master equation ODE, and the first passage time probability density (vesicle release time density) is directly and exactly computed from this ODE system as the transition rate to the final absorbing state. Therefore, Gillespie SSA simulations are not required, and the only Monte Carlo step involves drawing the vesicle release time from this exact probability distribution, without having to simulate/resolve any intermediate reaction times (once vesicle is released, one may have to recompute the FPT density, but that's not expensive). This should not affect the results, but would greatly improve the simulation efficiency, accuracy, and simplify the parameter sensitivity analysis. Since this should not affect the results, I would not request any significant modifications, but this could be somehow reflected in Materials and methods and checked (unless I misunderstand something and the transition rates are in fact stochastic).
1.3) Subsection “Rate equations of the simulated models”: this part of Materials and methods is particularly hard to read, and could be improved just by a simple restructuring. It is awkward to explain parameters before the actual equations are shown. The rate dependence on [Ca^{2+}] in subsection “Rate equations of the simulated models” seems strange: shouldn't there be a power of "n" in both terms in the denominator?
1.4) Discussion paragraph ten: when discussing alternative STF scenarios beyond "pVrbased" models, it could be appropriate to specifically distinguish between models where all vesicles have the same properties, and models with distinct vesicle pools with different properties (such as the superprimed pool models cited in various parts of the manuscript). Of course, this distinction is complicated, since vesicle pool heterogeneity could result from vesicles being in different transient states along the same slow priming processes, but this could also be pointed out… If the authors find it appropriate, they could also mention here the putative "highly Ca^{2+}sensitive pools" at endocrine cells, which I find interesting and potentially relevant (reviewed in Pedersen and Sherman (2009) PNAS 106:74327436).
1.5) It would be valuable to briefly comment on the facilitation decay time constant predicted by their model, since facilitation decay provides a powerful additional constraint on the potential model mechanism of STF. Even if experimental data on STF decay is not available in this case, the decay time constant could be presented as a prediction to be verified in future experimental work. The authors could simply quote the relevant priming/unpriming time scale determining STF decay to avoid additional simulation. I note that matching facilitation decay time course puts a strong constraint on the singlesensor and dualsensor models, and is an extra issue with such models among those mentioned in this work (see e.g. Matveev et al., 2002, cited in this manuscript).
1.6) Discussion paragraph nine: since the site activation model is not described in detail in Results or Discussion, but is one of two presented new models reproducing experimental observations, it would be helpful to describe it slightly more in this part of the manuscript, if not earlier (12 more sentences would suffice).
1.7) When the dualsensor model is first introduced, it would be appropriate to more prominently cite the related work of Sun et al., 2007, (which the authors cite elsewhere). Further, when discussing various STF models and prior work in this direction, I suggest citing the study of Ma et al., 2015, which is one of the most detailed models of STF and involves a fully stochastic approach: https://www.ncbi.nlm.nih.gov/pubmed/25210157.
2.1) While the characterization of the vesicle distribution is well done. It remains possible that the vesicles are not drawn independently from the distribution. In an extreme case, the number of vesicles could be fixed, or they could repulse each other, and yet the same distribution could still be found. I wonder if this can be discussed.
2.2) The role of heterogeneity in release and averaging across experiments (alluded to in Discussion paragraph six) is not clear to me. It would be good to know that the authors have convinced themselves that such heterogeneity does not overestimate variance or otherwise change the results.
3.1) The EM analysis was conducted in mutant animals. Exact genotypes were not provided, but based on the information in ReddyAlla et al., 2017 where the data were originally published I think they are working with Unc13 null animals expressing UASUnc13A via a Gal4 driver. This data is being used to drive models of normal synaptic function, so it's critical that this experiment is conducted in wildtype animals with normal levels of Unc13.
3.2) The authors use the data presented in Figure 1D to place synaptic vesicles relative to Ca^{2+} channels in their models. It's difficult to discern the number of animals, NMJs and active zones from the numbers provided: "n=19 observations in 10 EM cross sections/cells." This information should be clearly stated. Assuming an observation is a docked vesicle and a cross section is an active zone, the 3D placement of vesicles was derived from the observation of 19 docked synaptic vesicles in single sections of 10 active zones. Given that the central thesis that vesicles are heterogeneously distributed, so few observations cannot confidently capture the biological range.
This also suggests less than two docked vesicles per cross section, which is somewhat lower than docked vesicle/active zone numbers at wildtype Drosophila type Ib synapses previously reported by multiple groups, including the authors' prior work. This could be the genotype (see below) and/or the fact that single sections do not accurately capture the full complement of docked vesicles at an active zone. For example, vesicles that appear close to, but not in contact, with the membrane in one slice, may be in clear contact with the membrane in an adjacent plane. 3D EM approaches, which have been done at the Drosophila NMJ, would provide much better estimates of synaptic vesicle topology and obviate the need for deriving 3D estimates based on limited information.
3.3) It is well documented that the hundreds of active zones at Drosophila NMJs are of different developmental stages with very different morphologies and release properties, so the distribution of synaptic vesicles observed at a handful of active zones can't accurately represent the NMJ as a whole. This will be very challenging to address experimentally, but should at least be considered in their interpretations and addressed in the Discussion.
3.4) The study seems to involve a mixture of analysis of two motorneuron subtypes with different structural and functional properties without considered this in the modeling. Though not always stated, it looks like the EM, STED and GCaMP experiments were conducted at type Ib synapses, while the electrophysiology measured the compound response to both Ib and Is motor inputs. Can they specifically measure and model type Ib?
3.5) I think wildtype animals were investigated in the STED and electrophysiology experiments, but the genotypes are not noted. All genotypes should be clearly labeled in figures. Additionally, since the data is being reused here, this manuscript should provide all relevant methods rather than referring readers to the earlier publication. The STED results appear to be based on three animals from a single previously published experiment. Have any technical replicates of this experiment been performed to control for experimental variability?
3.6) There are very few references to other Drosophila labs working in this area. Multiple labs have conducted relevant work on the distribution of synaptic vesicles, Ca^{2+} channels, and heterogeneous release properties at this synapse that should be cited.
https://doi.org/10.7554/eLife.51032.sa1Author response
Essential revisions:
1.1) Equation 2: The whole Ca^{2+} current calibration is very hard to follow, and is my main concern regarding methods description. In Equation 2, shouldn't K_{M} be raised to power "h"? Also, this expression gives the relationship between fluorescence and intracellular Ca^{2+}, or was the fluorescence recorded with extracellular GCaMP? Further, why is K_{M} determined from GCaMP6 calibration? Isn't K_{M} an innate property of the Ca channel conductivity (and possibly Cadependent inactivation), having nothing to do with the indicator dye? The entire Ca^{2+} current calibration has to be explained much more clearly, I could not follow it.
The reviewer is correct that the parameter K_{M} should also be raised to the power of h. We have corrected this. Thank you for pointing this out. The relevant Ca^{2+} concentration is that of the extracellular medium, we also did not state this clearly in the equation.
We are sorry that it was unclear from the text what the reasoning of the experiment was. This is because we introduced an analysis of how Ca^{2+} influx depends on the extracellular Ca^{2+} concentration before explaining why this at all would be relevant. We have now restructured the text and swapped the sequence of Figures 2 and 3. This way we first present the electrophysiological recordings (to which the mathematical models are later contrasted) at different extracellular Ca^{2+} concentrations. This then delivers the explanation of why the dependency of Ca^{2+} influx a as a function of extracellular Ca^{2+} concentrations should be studied in Figure 3.
We now also extended the main text to clearly explain the premise of this experiment, which was (as the reviewer points out) to quantify the saturation behavior of the Ca^{2+} current as a function of the extracellular Ca^{2+} concentration. The reviewer is entirely correct that this is an intrinsic property of the Ca^{2+} channels. As such, it would be more direct to conduct an electrophysiology experiment to measure these currents, but this is not possible in this preparation. We therefore used an intracellular Ca^{2+} fluorescence measurement as a proxy. GCaMP6m was expressed within the presynapse and its relative fluorescence change in response to action potential activation (10 AP at 20 Hz) z) was measured. This only allows us to measure relative changes in Ca^{2+} inflow. We now point out more clearly that this experiment cannot derive absolute Ca^{2+} concentration because the system cannot be calibrated in such a way. The fluorescence values detect residual Ca^{2+} after the AP train. The relative change in fluorescence can give an indication of the relative magnitude of Ca^{2+} influx as a function of extracellular Ca^{2+} levels. We indeed see that the fluorescence increase saturates (as expected from the saturation of Ca^{2+} currents) and this is well described by a hill curve. We were most interested in the “IC50/KM” value of the curve, because we used this value to approximate the saturation of the Ca^{2+} current in Equation 4. The second relevant parameter, “h”, describes the cooperativity of the fluorescence increase as a function of extracellular Ca^{2+}. Because we reasoned that this should be a property of the GCaMP sensor, rather than the channel (this was described in Barnett et al., 2017) this parameter was fixed to the reported value of h = 2.43. The best fit parameter value of KM in this case was 2.679 mM, which is near identical to the KM value (2.6 mM) previously reported for the Ca^{2+} influx at the mammalian calyx of Held synapse (Schneggenburger and Neher, 1999). Both the IC50 and the Michaelis Menten constant report on the halfsaturation point. So as a matter of fact, our results are primarily confirmatory in that the Ca^{2+} channels at the Drosophila NMJ show a largely similar saturation behavior.
We also checked whether fixing the parameter “h” had a major effect on the estimation of KM from our experiments. This was not the case, as allowing additional variation of h during the fit resulted only in a marginally different estimate of the parameters (h=1.887 and K_{M}=3.054 mM). With the added explanation, we are now confident that the aim and rationale and the methodology is clear.
1.2) Unless I misunderstood the modeling/methods, the described stochastic simulation process, while technically correct, appears needlessly complicated and computationally expensive. Any Markov Chain with a single absorbing state and deterministic propensities ([Ca] is deterministic here) is described exactly by its master equation ODE, and the first passage time probability density (vesicle release time density) is directly and exactly computed from this ODE system as the transition rate to the final absorbing state. Therefore, Gillespie SSA simulations are not required, and the only Monte Carlo step involves drawing the vesicle release time from this exact probability distribution, without having to simulate/resolve any intermediate reaction times (once vesicle is released, one may have to recompute the FPT density, but that's not expensive). This should not affect the results, but would greatly improve the simulation efficiency, accuracy, and simplify the parameter sensitivity analysis. Since this should not affect the results, I would not request any significant modifications, but this could be somehow reflected in Materials and methods and checked (unless I misunderstand something and the transition rates are in fact stochastic).
In principle this is correct, however, to calculate the solution of the master equation ODE can only be done numerically, due to the time varying Ca^{2+}dependent transition rates. Then the FPT density has to be calculated from this numerical distribution. Moreover, this distribution would have to be calculated repeatedly, for each randomly drawn SV:channel distance, and after each fusion, where first a random time for replenishment has to be drawn, then the numerical solution of the ODE has to be found to calculate the FPT density. We doubt that this will improve simulation efficiency or accuracy.
1.3) Subsection “Rate equations of the simulated models”: this part of Materials and methods is particularly hard to read, and could be improved just by a simple restructuring. It is awkward to explain parameters before the actual equations are shown. The rate dependence on [Ca^{2+}] in subsection “Rate equations of the simulated models” seems strange: shouldn't there be a power of "n" in both terms in the denominator?
We agree on this point and the section has been restructured. The equations were moved up before the parameter explanations. We added: “Parameters are described below. For designation and value of parameters, see Table 2 and 3.”
The reviewer is of course right about the power of ‘n’ on the terms in the denominator. We doublechecked that this was correctly implemented in the scripts used for simulations and the typo in the manuscript has been corrected.
1.4) Discussion paragraph ten: when discussing alternative STF scenarios beyond "pVrbased" models, it could be appropriate to specifically distinguish between models where all vesicles have the same properties, and models with distinct vesicle pools with different properties (such as the superprimed pool models cited in various parts of the manuscript). Of course, this distinction is complicated, since vesicle pool heterogeneity could result from vesicles being in different transient states along the same slow priming processes, but this could also be pointed out… If the authors find it appropriate, they could also mention here the putative "highly Ca^{2+}sensitive pools" at endocrine cells, which I find interesting and potentially relevant (reviewed in Pedersen and Sherman (2009) PNAS 106:74327436).
This is a good point, and we have inserted such a discussion now, distinguishing between vesicle pools with different intrinsic properties, and vesicles with the same intrinsic properties, but different distances to the Ca^{2+}channels. Regarding the highly Ca^{2+}sensitive pool, we are well aware of the work and have cited it repeatedly. However, it is based on work in chromaffin cells and beta cells, which is now 15 years old. We would therefore prefer not to go into it here, as we fear it might sidetrack the text.
1.5) It would be valuable to briefly comment on the facilitation decay time constant predicted by their model, since facilitation decay provides a powerful additional constraint on the potential model mechanism of STF. Even if experimental data on STF decay is not available in this case, the decay time constant could be presented as a prediction to be verified in future experimental work. The authors could simply quote the relevant priming/unpriming time scale determining STF decay to avoid additional simulation. I note that matching facilitation decay time course puts a strong constraint on the singlesensor and dualsensor models, and is an extra issue with such models among those mentioned in this work (see e.g. Matveev et al., 2002, cited in this manuscript).
This is a very good suggestion. We now mention this reference in this context in the text. We have also followed the reviewer’s suggestion to predict the timecourse of facilitation of this model which is now included as a supplementary item (Figure 7—figure supplement 3). These predictions can then later be tested experimentally by varying the ISI in PPR experiments.
1.6) Discussion paragraph nine: since the site activation model is not described in detail in Results or Discussion, but is one of two presented new models reproducing experimental observations, it would be helpful to describe it slightly more in this part of the manuscript, if not earlier (12 more sentences would suffice).
Thank you for pointing this out. We have now added a more elaborate description of the site activation model in the Results section at the end of the paragraph “Rapidly regulating the number of participating release sites accounts for eEJC_{1} amplitudes, STF, temporal transmission profiles and variances”.
1.7) When the dualsensor model is first introduced, it would be appropriate to more prominently cite the related work of Sun et al., 2007, (which the authors cite elsewhere). Further, when discussing various STF models and prior work in this direction, I suggest citing the study of Ma et al., 2015, which is one of the most detailed models of STF and involves a fully stochastic approach: https://www.ncbi.nlm.nih.gov/pubmed/25210157.
Thanks for pointing this out. We have added statements and references accordingly.
2.1) While the characterization of the vesicle distribution is well done. It remains possible that the vesicles are not drawn independently from the distribution. In an extreme case, the number of vesicles could be fixed, or they could repulse each other, and yet the same distribution could still be found. I wonder if this can be discussed.
There are two issues here: one is what the vesicle distribution looks like, the other is how this distribution comes about. The first question we can solve (and we believe we have solved it), the second question is open. However, for the work we present the question how the distribution comes about (by independent random localization of each vesicle, or affected by repulsive interaction between vesicles) does not affect the initial positioning of vesicles – the random drawing we employ will ensure a distribution of vesicles, which closely reflects the identified distribution, even if the mechanism that created this distribution does not involve random positioning. The only place where the mechanism creating the distribution becomes a question is during recovery, i.e. when some vesicles have fused as a result of the first AP and new vesicles need to prime. It is possible that vesicle positions are again drawn independently from the distribution, alternatively vesicles could selectively prime into the now empty positions (‘slots’). We have chosen the latter possibility for several reasons: (1) It is the only mechanism that would have any chance of producing facilitation in pVrbased models. The vesicles fusing during the first AP are localized much closer to the Ca^{2+}channels than the average vesicle, and priming into newly selected random positions will therefore shift the distribution of vesicles away from the Ca^{2+}channels, making facilitation an impossibility. As we conclude that facilitation using pVrbased models is not likely, we had to make sure we were not stacking the deck against this possibility; (2) If vesicles prime into new slots, the distribution of vesicles after full recovery would have shifted after the action potential, and then would have to shift back, which is a complicated mechanism for which there is no evidence; (3) We believe it is generally accepted that in synapses vesicles prime into preexisting ‘slots’ on the plasma membrane, created by priming proteins (including Unc13s). And we have shown that these proteins are extremely stable (in FRAP experiments we previously established a recovery time constant of ~6 hours, ReddyAlla et al., 2017). Therefore, it is effectively the slots we have to draw according the distribution, and not the vesicle positions, although the position of vesicles of course follows from the distribution of slots.
2.2) The role of heterogeneity in release and averaging across experiments (alluded to in Discussion paragraph six) is not clear to me. It would be good to know that the authors have convinced themselves that such heterogeneity does not overestimate variance or otherwise change the results.
We have now made the text clearer. In this kind of experiment there are two sources of variance: there is the variance between responses within the same cell (withincell variance), and there is variance between different cells (betweencell variance). An advantage of having a stochastic model is that we can compare it not only to the mean response, but also to the variance. The relevant variance is the withincell variance, as it is measured experimentally in variancemean experiments, where a single cell is stimulated repetitively. Therefore, we have displayed the data as the mean (between cells) plusminus the mean withincell variance in the figures where simulation output is compared to experiments (Figures 4, 6, 7). In fact, this results in a very good match between model and data, both for the means and the variances, for the depriming model (e.g. Figure 7E, G). This does not overestimate variance. Indeed, a single model should be able to reproduce the withincell variance, but it cannot encompass the betweencell variance, since we would need one model per cell.
We also tried to clarify the statement referred to by replacing the second part of the statement with the following sentence: “Note that such parabola agrees with the meanvariance relationship in a binomial distribution. However these simulations show that the assumption of heterogeneous release probability can also lead to the experimentally observed parabolic variancemean relationship.”
3.1) The EM analysis was conducted in mutant animals. Exact genotypes were not provided, but based on the information in ReddyAlla et al., 2017 where the data were originally published I think they are working with Unc13 null animals expressing UASUnc13A via a Gal4 driver. This data is being used to drive models of normal synaptic function, so it's critical that this experiment is conducted in wildtype animals with normal levels of Unc13.
We are sorry that this was not clearly stated and now provide all details in a section of the Materials and methods where all genotypes are listed by figures so this information can easily be found. The reviewer is correct that this is a previously published dataset and the flies are indeed Unc13 null animals expressing the Unc13A isoform under an elavGal4 neuronal driver. The rationale here was to isolate the contribution of Unc13A, because the second Drosophila isoform Unc13B only marginally (<5%) contributes to synaptic transmission. As Unc13B localizes at larger distances from the Ca^{2+} channels and likely functions at docking at those distant positions (which is hard to map in the EM images) we thought that this would be more adequate. Also, the second approach to map the docking sites using STED microscopy only maps the distribution of the Un13A isoform (it is an isoform specific antibody targeting the Nterminal region of the protein). The reviewer is correct that the reexpression of the protein using the elavGal4 driver could potentially lead to a higher expression than the endogenous promoter. However, we would like to point out, that our focus here is to study the distribution of docking sites rather than their number (see also point 3.2). Therefore, the higher levels of Unc13A would not affect the conclusions of our study as long as the protein shows a similar distribution. Nonetheless, the reviewer makes a valid point that a corresponding analysis should also be performed in wildtype animals. We therefore now added this analysis to our study which is also based on a previously published dataset (measurements of docked SVs in wildtype Drosophila larvae from Böhme et al., 2016.) The result and the fitted Rayleigh distribution are shown in Figure 1—figure supplement 1A.
When comparing the relevant (2D) integrated Rayleigh distribution from which vesicles are placed in our simulation in both cases (wildtype vs. Unc13A rescue, Figure 2—figure supplement 2B) it becomes evident that the distributions are near identical. This rules out that the altered Unc13A levels severely affect the docked vesicle distribution.
This alleviates the reviewer’s concern that the genotype influenced our analysis of docked SV positions. Please note that we did not rerun the optimization of all models with the second distribution. This would be pointless, because the distributions are essentially identical and any effect on the best fit parameters would be marginal and it would not be justified to spend weeks of calculations on a computer grid for this purpose.
3.2) The authors use the data presented in Figure 1D to place synaptic vesicles relative to Ca^{2+} channels in their models. It's difficult to discern the number of animals, NMJs and active zones from the numbers provided: "n=19 observations in 10 EM cross sections/cells." This information should be clearly stated. Assuming an observation is a docked vesicle and a cross section is an active zone, the 3D placement of vesicles was derived from the observation of 19 docked synaptic vesicles in single sections of 10 active zones. Given that the central thesis that vesicles are heterogeneously distributed, so few observations cannot confidently capture the biological range.
We are sorry that this information was not provided here. We indeed only referred to the original publication. Of course it makes sense to immediately make the information available and we have now added this to the figure legend. We also provide all source data and provide the relevant information in the transparent reporting form. The analysis for the Unc13A reexpression dataset is based on 10 AZs from at least two animals. The analysis of the EM data from wildtype flies (Figure 1—figure supplement 1A) is based on 11 AZs from 5 animals.
We disagree with the reviewer regarding the second point (that the number of observations is too low to conclude on the distribution of distances) and the argument against this is provided below. In this study, we find that the distribution of distances is broad and this can be safely concluded from the number of observations we have. We used two different methods to test how reliable the finding of a broad vesicle distance distribution from the EM data is: by using confidence intervals from the fitting of the Rayleigh distribution presented in the paper (based on the likelihood function) and by bootstrapping from the data set at hand (without assuming anything about the underlying distribution).
Likelihood approach:
In the paper, we fitted a Rayleigh distribution to the data set using the maximum likelihood function, which yielded σ=76.52 nm as the best fit parameter.
From the likelihood function we can estimate 95% confidence intervals of σ and the means and SDs of the found distributions (Author response table 1).
Using the three different σvalues from above (best fit and the lower and upper bound of the 95% CI), the integrated Rayleigh distribution from which we draw SV distances in our simulations looks as shown in Author response image 1:
Since the STD varies with the choice of parameter in a monotone way and since all these three distributions are relatively broad, any choice of parameter in the confidence interval will yield a broad distribution. Thus, when fitting a Rayleigh distribution to the data, we get a broad distribution of distances, also when taking the uncertainty of the fit into account.
Bootstrap analysis from the data
To estimate the mean and SD of the underlying distribution of SV distances without making any assumptions on the nature of this distribution we bootstrapped from the data set at hand (Figure 1B, n=19). We drew 19 observations with replacement from the data set 10000 times and estimated the mean and the variance in each set of observations. The estimated mean is the average sampled mean, and the estimated SD is the square root of the average sampled variance. The confidence intervals are the 2.5% and 97.5% quantiles. In the case of the SD, the confidence interval is the square root of the quantiles of the variances.
The histogram in Author response image 2 shows the sampled mean with 95% confidence interval marked:
In this way, we use the knowledge we have about the data (our data set of 19 observations) to estimate the underlying mean and SD without making any assumptions about the distribution. Author response table 2 compares the mean and SD from bootstrapping to the mean and SD found by fitting the Rayleigh distribution (i.e. not the integrated distribution, since we only estimate the mean and SD of the distribution underlying our 1D sampling):
Thus, both methods of estimation indicate that the underlying distribution is in fact broad. The first method provides a confidence interval for the distribution fitted, whereas the bootstrap estimates mean and SD (and confidence intervals) from the number of observations in the dataset without assuming the underlying distribution to be a Rayleigh distribution (i.e. it is entirely nonparametric). Since both methods yield a broad distribution (large SD) and since the integration of the distribution will lead to a larger SD, both methods points toward the SV distance distribution being broad.
As the reviewer pointed out, the broadness of the distribution is at the core of the problem investigated in this paper, and we think this analysis justifies our assumptions on this matter. Besides, the distribution is confirmed in an independent EM dataset (Figure 1—figure supplement 1A) and in independent STED analyses (Figure 1E,F, Figure 1—figure supplement 1C,D).
This also suggests less than two docked vesicles per cross section, which is somewhat lower than docked vesicle/active zone numbers at wildtype Drosophila type Ib synapses previously reported by multiple groups, including the authors' prior work. This could be the genotype (see below) and/or the fact that single sections do not accurately capture the full complement of docked vesicles at an active zone. For example, vesicles that appear close to, but not in contact, with the membrane in one slice, may be in clear contact with the membrane in an adjacent plane. 3D EM approaches, which have been done at the Drosophila NMJ, would provide much better estimates of synaptic vesicle topology and obviate the need for deriving 3D estimates based on limited information.
We realize that the image previously shown in Figure 1A was maybe not representative and included a new one. We would like to point out that we are not really interested in the absolute number of docked vesicles in this study (in fact, this is a free parameter in the model). The focus was primarily on their distribution. Nevertheless, the number of docked SVs does in fact align with previous estimates by us and others: the Unc13A reexpression dataset reports ~2 docked SV per active zone, the wildtype dataset ~3. Estimates by other groups are in this range, see e.g. Bruckner et al., 2016 (http://doi.org/10.1083/jcb.201601098) where the O’ConnerGiles lab estimated ~2.5 docked SVs/AZ in thin EM sections and 23 docked SVs/AZ from EM tomography images. One should note, however, that the definition of docking often depends on docking definitions and on the resolution in the EM images, and as such is not always directly comparable. However, as pointed out above, since this number is a free parameter for our modeling; it is not a concern that some SVs docked at this AZ may be missed (that would be docked in other sections).
Using EMtomography could be another method to address the docked SV distribution, but as such an approach relies heavily on computational approaches for image reconstruction, we feel that this is less direct. EMtomography might make it possible to identify more docked vesicles in a given presynapse; however, this will not change the distribution of docking sites. Besides, we show a complementary approach by our STED analysis, which we feel is stronger because it uses another technique but delivers a similar estimate of SV docking sites with molecular information.
3.3) It is well documented that the hundreds of active zones at Drosophila NMJs are of different developmental stages with very different morphologies and release properties, so the distribution of synaptic vesicles observed at a handful of active zones can't accurately represent the NMJ as a whole. This will be very challenging to address experimentally, but should at least be considered in their interpretations and addressed in the Discussion.
We fully agree that the heterogeneity across AZs is a very important factor and indeed inspired us to this project. Therefore it seems that this point may partly rely on a misunderstanding, as it is exactly the main premise of our study to quantitatively investigate the consequence of heterogeneous SV placemen. Note that in our model we make no distinction between whether the variability is present within single AZs, or is a variability between AZs, therefore the model captures both. Thus, our approach is actually suitable to address exactly the type of heterogeneity between AZs that the reviewer is pointing to. We feel our approach is much more accurate than previous ones by taking into account this heterogeneity and by furthermore doing so in quantitative manner (e.g. captured by a single parameter of the Rayleigh distribution) and showing its consequence for function. It is correct that apart from to this there are additional sources of heterogeneity and we now extended on the discussion of these sources (e.g. Ib/Is type boutons, heterogeneous priming).
3.4) The study seems to involve a mixture of analysis of two motorneuron subtypes with different structural and functional properties without considered this in the modeling. Though not always stated, it looks like the EM, STED and GCaMP experiments were conducted at type Ib synapses, while the electrophysiology measured the compound response to both Ib and Is motor inputs. Can they specifically measure and model type Ib?
We agree that many things can be learned here regarding additional sources of heterogeneity (see also point 3.3). We followed the reviewer’s suggestion and added this to in the Discussion. Indeed, both connection types likely contribute to the electrophysiological measurements (though it is currently not exactly known to which extent). Also, the STED analysis focusses on AZs within Ib type boutons and it could be interesting to systematically investigate the different bouton types using STED microscopy in the future. The GCaMP data were obtained from type Ib synapses, but measurements were merely used to capture the saturation behavior of Ca^{2+} currents as a function of extracellular Ca^{2+} concentrations (see also point 1.1). We found exactly the same relationship as described for the murine calyx of Held synapse. This makes sense, because this should be a fundamental property of the Ca^{2+} channels. We therefore feel that it is less likely this property would be markedly different between the two input types as it was so similar across species. Nevertheless, this could be investigated in the future. EM analysis favors Ib boutons (AZs are typically used from large boutons with prominent postsynaptic SSR). However, it cannot be excluded that sometimes Is connections are also sampled (because they are difficult to distinguish at this level).
We agree with the reviewer that mapping possible differences in function, composition and SV distribution might be interesting. We also agree that this will be very challenging, likely requiring the genetic manipulation of one vs. the other. This again could perturb the system in unexpected ways in addition to desired effects, so specificity could be a problem. Depending on what comparative analyses reveal, this could lead to a refinement of the model to the two bouton types (and how each of these additionally contributes to heterogeneity). However, this is beyond the scope of the current study. As it stands, the current model is a compound model describing overall synaptic transmission to this postsynaptic connection, which we want to publish as such.
3.5) I think wildtype animals were investigated in the STED and electrophysiology experiments, but the genotypes are not noted. All genotypes should be clearly labeled in figures. Additionally, since the data is being reused here, this manuscript should provide all relevant methods rather than referring readers to the earlier publication. The STED results appear to be based on three animals from a single previously published experiment. Have any technical replicates of this experiment been performed to control for experimental variability?
We now provide all details on genotypes in a dedicated section of the Materials and methods where genotypes are listed by figures, making this information immediately accessible. Electrophysiological recordings were performed in animals expressing an Ok6Gal4 driver: Ok6GAL4/+ (Ok6Gal4/II crossed to w1118). The reason that Gal4 driver was introduced into the wildtype strain was such that these animals could function as a control group for a motor neuron specific knockdown of Syt7 (which was suggested to function as a facilitation sensor). However, in the end we did not include these data (where the Ok6Gal4 is combined with a UASSyt7RNAi).
Technically, there could be a potential effect of the Ok6Gal4 driver on synaptic transmission, its Ca^{2+} dependence and/or its shortterm plasticity. To alleviate this concern, we performed electrophysiological recordings in Drosophila larva expressing Ok6Gal4 and compared these to wildtype animals (+/+: w1118) in parallel. These data revealed no difference between the two groups and are now shown in the new Figure 2—figure supplement 3.
We found no difference in the amplitude of APevoked eEJCs at 0.75 mM [Ca^{2+}]_{ext}: wildtype eEJC amplitude: 35.13 ± 24.94 nA (n = 6 cells from 3 animals, mean ± SD) and Ok6Gal4/+ eEJC amplitude: 32.36 ± 9 nA, (n = 8 cells from 4 animals, mean ± SD), P = 0.77; Student´s t – test. Similar responses were also seen at 1.5 mM extracellular Ca^{2+}, arguing for a similar Ca^{2+} sensitivity: wildtype eEJC amplitude 98.24 ± 10.73 nA (mean ± SD, n = 5 cells from 3 animals) and Ok6Gal4/+ eEJC amplitudes: 93.22 ± 12.16 nA (mean ± SD, n = 5 cells from 3 animals); P = 0.55; Student´s t – test.
Finally, both groups showed indistinguishable PPR values at both concentrations (0.75 mM Ca^{2+}: wildtype: 1.02 ± 0.51 (n = 6 cells from 3 animals, mean ± SD) and Ok6Gal4/+: 1.15 ± 0.26 (n = 8 cells from 4 animals, mean ± SD); P = 0.55; Student´s t – test.; 1.5 mM Ca^{2+}: wildtype: 0.65 ± 0.12 (mean ± SD, n = 5 cells from 3 animals) and Ok6Gal4/+: 0.58 ± 0.08 (mean ± SD, n = 5 cells from 3 animals); P = 0.39; Student´s t – test):
Together, these observations rule out that the presence of the Ok6Gal4 has a severe effect on the magnitude or Ca^{2+} dependence of APevoked synaptic transmission, or that it markedly alters the synapse’ shortterm plasticity characteristics in a way that could affect the conclusions of our study.
We furthermore included a description of the full procedure of larval preparation, fixation, immunohistochemistry, image acquisition and the averaging of the STED images in the Materials and methods section. And we also now include another biological replicate of the approach (completely different animals, different staining, imaged separately) based on 583 individual AZ images from 3 wildtype larva for comparison. This analysis essentially reproduces the distribution Unc13A distribution found before. These data are depicted a new figure, Figure 1—figure supplement 1.
3.6) There are very few references to other Drosophila labs working in this area. Multiple labs have conducted relevant work on the distribution of synaptic vesicles, Ca^{2+} channels, and heterogeneous release properties at this synapse that should be cited.
We apologize that this was missed and now provide additional references to these important works.
https://doi.org/10.7554/eLife.51032.sa2Article and author information
Author details
Funding
Deutsche Forschungsgemeinschaft (Emmy Noether Programme, Project Number 261020751)
 Alexander M Walter
Deutsche Forschungsgemeinschaft (Project Number 278001972  TRR 186)
 Alexander M Walter
Independent Research Fund Denmark (Pregraduate scholarship (814100007B))
 Jakob Balslev Sørensen
 Janus RL Kobbersmed
NeuroCure Cluster of Excellence (NeuroCure Fellowship)
 Andreas T Grasskamp
Einstein Stiftung Berlin (Einstein Center for Neuroscience)
 Meida Jusyte
 Alexander M Walter
University of Copenhagen (Data Science Laboratory)
 Janus RL Kobbersmed
Lundbeck Foundation (R2772018802)
 Jakob Balslev Sørensen
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
JRL Kobbersmed was supported by a pregraduate Scholarship in the Medical Sciences from the Independent Research Fund Denmark (application by JB Sørensen) and by the Data Science Laboratory, Faculty of Science, University of Copenhagen. AT Grasskamp was supported by a NeuroCure Ph.D. fellowship funded by the NeuroCure cluster of excellence within the International Graduate Program Medical Neurosciences (Charité Universitätsmedizin Berlin, Germany). AM Walter is a member of the Einstein Center for Neurosciencies Berlin, M Jusyte was supported by an Einstein Center for Neurosciences Ph.D. fellowship. This work was supported by grants from the Deutsche Forschungsgemeinschaft to AM Walter (Emmy Noether program project number 261020751, TRR 186 project number 278001972). We thank Matthijs Verhage and Niels Cornelisse for helpful comments on our manuscript. We would like to thank Victor Matveev for helpful comments on the use of his Ca^{2+} simulator software. We thank Barth van Rossum, LeibnizForschungsinstitut für Molekulare Pharmakologie (FMP), for the illustrations in Figure 9. We also thank The Bioinformatics Centre, University of Copenhagen, for the use of their servers for heavy computations. Stocks obtained from the Bloomington Drosophila Stock Center (NIH P40OD018537) and from the Kyoto Stock Center were used in this study.
Senior Editor
 Ronald L Calabrese, Emory University, United States
Reviewing Editor
 Mark CW van Rossum, University of Nottingham, United Kingdom
Reviewer
 Victor Matveev, NJIT
Publication history
 Received: August 12, 2019
 Accepted: February 14, 2020
 Accepted Manuscript published: February 20, 2020 (version 1)
 Version of Record published: April 9, 2020 (version 2)
Copyright
© 2020, Kobbersmed 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

 1,928
 Page views

 440
 Downloads

 1
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.