We first set out to quantify the SV release site:Ca²⁺ channel topology. For this we analysed EM micrographs of AZ cross-sections and quantified the distance between docked SVs (i.e. SVs touching the plasma membrane) and the centre of electron dense ‘T-bars’ (where the voltage gated Ca²⁺ 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 (‘EM-dataset Unc13A rescue’) (Reddy-Alla 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 cross-section of a three-dimensional synapse. To derive the relevant coupling distance distribution for all release sites (including the ones outside the cross-section), the Rayleigh distribution was integrated around a circle (Figure 1C), resulting in the following probability density function (pdf, see Materials and methods for derivation):

Figure 1 with 1 supplement see all

Download asset Open asset

Figure 1—source data 1

Source data for graphs in Figure 1 and Figure 1—figure supplement 1.

https://cdn.elifesciences.org/articles/51032/elife-51032-fig1-data1-v2.xlsx

Download elife-51032-fig1-data1-v2.xlsx

Figure 1—source data 2

Matlab codes used for data analysis, original images, and instructions for analysis depicted in Figure 1 and Figure 1—figure supplement 1.

https://cdn.elifesciences.org/articles/51032/elife-51032-fig1-data2-v2.zip

Download elife-51032-fig1-data2-v2.zip

These pdfs were more symmetrical than the ones from the cross-sections 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²⁺ channel coupling distances without relying on the integration of docked SV observations from cross-sections: since (M)Unc13 was recently described as a molecular marker of SV release sites (Reddy-Alla 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 ring-like distribution of the Unc13A fluorescence. In previous works we had established that the voltage gated Ca²⁺ 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).

Having identified the high degree of heterogeneity in the docked SV:Ca²⁺ channel coupling distances, we became interested in how this affected synaptic function. We therefore characterized synaptic transmission at control NMJs (Ok6-GAL4 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²⁺ concentration of the extracellular solution which affects AP-induced Ca²⁺ influx (see below). We used this approach and investigated responses evoked by repetitive (paired-pulse) 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₁ amplitudes) with increasing extracellular Ca²⁺ (Figure 2A,B). STP was assessed by determining the paired-pulse ratio (PPR): the amplitude of the second response divided by first. The eEJC₂-amplitude was determined taking the decay of eEJC₁ into account (see insert in Figure 2C, Figure 2—figure supplement 1A). At low extracellular Ca²⁺ (0.75 mM), we observed strong STF (with an average PPR value of 1.80), which shifted towards depression (PPR < 1) with increasing Ca²⁺ concentrations (Figure 2C,D). Thus, the same NMJ displays both facilitation and depression depending on the extracellular Ca²⁺ concentration, making this a suitable model synapse to investigate STP behaviour.

Figure 2 with 3 supplements see all

Download asset Open asset

Figure 2—source data 1

Raw data for experiments displayed in Figure 2 and Figure 2—figure supplement 2.

https://cdn.elifesciences.org/articles/51032/elife-51032-fig2-data1-v2.zip

Download elife-51032-fig2-data1-v2.zip

Figure 2—source data 2

Matlab code used for data analysis of electrophysiological traces.

https://cdn.elifesciences.org/articles/51032/elife-51032-fig2-data2-v2.zip

Download elife-51032-fig2-data2-v2.zip

Figure 2—source data 3

Raw data which was used for depicted anaylsis in Figure 2 and Figure 2—figure supplement 3.

https://cdn.elifesciences.org/articles/51032/elife-51032-fig2-data3-v2.xlsx

Download elife-51032-fig2-data3-v2.xlsx

In panels B and D the mean eEJC₁ 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 AP-evoked responses at the same NMJ between trials (10 s apart) at different extracellular Ca²⁺ concentrations (Figure 2E,F). At low concentrations (0.75 mM), the probability of transmitter release is low, resulting in a low mean eEJC₁ amplitude with little variation (Figure 2E,F, Figure 2—figure supplement 2 ). With increasing extracellular Ca²⁺, the likelihood of SV fusion increased and initially so did the variance (e.g. at 1.5 mM extracellular Ca²⁺). However, further increase in extracellular Ca²⁺ (3 mM, 6 mM, 10 mM) led to a drop in variance (Figure 2E, Figure 2—figure supplement 2). Figure 2F depicts this average ‘variance-mean’ 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).

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²⁺ concentration (Neher and Sakaba, 2008). Secondly, because of the strong buffering of Ca²⁺ signals at the synapse, the magnitude of the AP-evoked Ca²⁺ transients dramatically declines with distance from the Ca²⁺ 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²⁺ 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 AP-induced fusion events of docked SVs placed according to the found distribution. A prerequisite for this is to first faithfully simulate local, AP-induced Ca²⁺ 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²⁺ channels in the base centre (Meinrenken et al., 2002). The base dimensions (length = width) were determined as the mean inter-AZ distance of all AZs to their four closest neighbours (because of the 4-fold symmetry) from NMJs stained against the AZ-marker 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²⁺ channel is the only relevant parameter) covering the same AZ area (Figure 3B, Table 1).

Figure 3 with 2 supplements see all

Download asset Open asset

Figure 3—source data 1

Average cell-wise mean fluorescence values and fit parameters of hill-curve fit on presynaptic GCaMP data, and NND values.

https://cdn.elifesciences.org/articles/51032/elife-51032-fig3-data1-v2.xlsx

Download elife-51032-fig3-data1-v2.xlsx

To simulate the electrophysiological experiments above, where the extracellular Ca²⁺ concentration was varied (Figure 2), it was important to establish how the extracellular Ca²⁺ concentration influenced AP-induced Ca²⁺ influx. In particular, it is known that Ca²⁺ currents saturate at high extracellular Ca²⁺ 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 fluorescence-based approach as a proxy. AP-evoked Ca²⁺ influx was assessed in flies presynaptically expressing the Ca²⁺-dependent fluorescence reporter GCaMP6m (;P{y[+t7.7] w[+mC]=20XUAS-IVS-GCaMP6m}attP40/Ok6-GAL4). Fluorescence increase was monitored upon stimulation with 20 APs (at 20 Hz) while varying the extracellular Ca²⁺ concentration and showed saturation behaviour for high concentrations (Figure 3—figure supplement 1). This is consistent with a previously described Michaelis-Menten type saturation of fluorescence responses of a Ca^2+-sensitive dye upon single AP stimulation at varying extracellular Ca²⁺ concentrations at the Calyx of Held, where half-maximal Ca²⁺ influx was observed at 2.6 mM extracellular Ca²⁺ (Schneggenburger et al., 1999). This relationship was successfully used in the past to predict Ca²⁺ 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²⁺ and therefore used this value as K_M,current in a Michaelis-Menten equation (Materials and methods, Equation 5) to calculate AP-induced presynaptic Ca²⁺ influx. The second parameter of the Michaelis-Menten equation, (the maximal Ca²⁺ 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²⁺ concentrations at rest were also slightly dependent on the extracellular Ca²⁺ levels in a Michaelis-Menten relationship with the same dependency (K_M,current) and a maximal resting Ca²⁺ concentration of 190 nM (resulting in 68 nM presynaptic basal Ca²⁺ concentration at 1.5 mM external Ca²⁺). With these and further parameters taken from the literature on Ca²⁺ diffusion and buffering (see Table 1) the temporal profile of Ca²⁺ 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.

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 well-behaved 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; Vere-Jones, 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 inter-animal 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²⁺ binding states based on inhomogeneous, continuous time Markov models with transition rates governing reaction probabilities (see Materials and methods for details).

Figure 4 with 1 supplement see all

Download asset Open asset

Figure 4—source data 1

Simulation data for graphs in Figure 4C, G, E and Figure 4—figure supplement 1.

https://cdn.elifesciences.org/articles/51032/elife-51032-fig4-data1-v2.xlsx

Download elife-51032-fig4-data1-v2.xlsx

We also needed to consider where new SVs would (re)dock once SVs had fused and implemented the simplest scenario of re-docking in the same positions. This ensures a stable distribution over time and agrees with the notion that vesicles prime into pre-defined release sites, which are stable over much longer time than a single priming/unpriming event (Reddy-Alla et al., 2017).

The first model we tested was the single-sensor model proposed by Lou et al. (2005), where an SV binds up to 5 Ca²⁺ 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²⁺ (un)binding kinetics were taken from Wölfel et al. (2007) Table 3, the values of the maximal Ca²⁺ 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₁ amplitudes, and the mean and variance of PPRs at various extracellular Ca²⁺ concentrations and contrasted these with the experimental data.

This single-sensor model was able to reproduce the eEJC₁ values (Figure 4B,C). Moreover, the model accounted for the STD typically observed at high extracellular Ca²⁺ 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²⁺ (Figure 4D,E) and the variances predicted by this model were much larger than found experimentally (Figure 4F,G). The observation that eEJC₁ 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²⁺-triggered release event (Lou et al., 2005). As this model lacks a specialized mechanism to induce facilitation, residual Ca²⁺ binding to the Ca²⁺ 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²⁺ 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 paired-pulse experiment at 0.75 mM extracellular Ca²⁺ (Figure 5).

Figure 5

Download asset Open asset

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₁, immediately after AP₁, immediately before AP₂ (i.e. after refilling) and immediately after AP₂ (the external Ca²⁺ concentration was 0.75 mM). From this analysis it becomes clear that the pV_r1 caused by AP₁ 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 high-pV_r SVs are depleted by AP₁. SV replenishment refills the majority (but not all) of those sites and thus AP₂/pV_r2 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²⁺, 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²⁺ 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.

The single-sensor model failed to reproduce the experimentally observed STF at low extracellular Ca²⁺ concentrations because of the dominating depletion of SVs close to Ca²⁺ channels, and the inability to draw on SVs further away. However, this situation may be improved by a second Ca²⁺ sensor optimized to enhance the pV_r2 in response to AP₂. Indeed, in the absence of the primary Ca²⁺ sensor for fusion, Ca²⁺ sensitivity of synaptic transmission persists, which was explained by a dual sensor model (Sun et al., 2007). It was recently suggested that syt-7 functions alongside syt-1 as a Ca²⁺ sensor for release (Jackman et al., 2016), and deterministic mathematical dual fusion-sensor model assuming homogeneous release probabilities (which implies homogeneous SV release site:Ca²⁺ 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 fusion-sensor model is that while syt-1 is optimized to detect the rapid, AP-induced Ca²⁺ transients (because of its fast Ca²⁺ (un)binding rates, but fairly low Ca²⁺ affinity), the cooperating Ca²⁺ sensor is optimized to sense the residual Ca²⁺ after this rapid transient (Figure 3C) (with slow Ca²⁺ (un)binding, but high Ca²⁺ affinity). The activation of this second sensor after (but not during) AP₁ could then enhance the release probability of the remaining SVs for AP₂ (Figure 6A,B). This is illustrated in Figure 6B, where k₂ (the on-rate of Ca²⁺ binding to the slow sensor) is varied resulting in different time courses and amounts of Ca²⁺ 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²⁺ ion (Figure 6A), resulting in multiplicative effects on the SV fusion rate. While the fast fusion reaction appears to have a 5-fold Ca²⁺ cooperativity (Bollmann et al., 2000; Burgalossi et al., 2010; Schneggenburger and Neher, 2000), it is less clear what the Ca²⁺ cooperativity of a second Ca²⁺ 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²⁺ 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 syt-1 /-2), or whether it is present at the plasma membrane. Both scenarios are functionally possible and it was indeed reported that syt-7 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²⁺ affinity of K_D = 1.5 μM (Brandt et al., 2012; Jackman and Regehr, 2017).

Figure 6 with 1 supplement see all

Download asset Open asset

Figure 6—source data 1

Simulation data for graphs in Figure 6B, E, G, I and Figure 6—figure supplement 1B, D, E.

https://cdn.elifesciences.org/articles/51032/elife-51032-fig6-data1-v2.xlsx

Download elife-51032-fig6-data1-v2.xlsx

Like for the single-sensor 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 single-sensor model) together with k₂ (Ca²⁺ 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₂ had an effect on the PPR in simulations, confirming that the second sensor was able to improve the release following AP₂ (Figure 6C). Figure 6D–I show that the dual fusion-sensor model could fit the eEJC₁ amplitudes and the model slightly improved the higher PPR values at the low- and the lower PPR values at high extracellular Ca²⁺ 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²⁺were ~ 1.08 in the simulation compared to ~ 1.80 in the experiments). Another problem of the dual fusion-sensor model was that release became more asynchronous than observed experimentally (Figure 6D), which was due to the triggering of SV fusion in-between 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²⁺ 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₁ amplitudes (Figure 6K). Changing s had the largest effect when k₂ 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₁ (when increasing k₂) or by decreasing the effect on AP₂ (when decreasing k₂), which both counteracts STF (Figure 6B,J).

Fitting the dual fusion-sensor model with a Ca²⁺ 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 fusion-sensor 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 syt-7 functions in STF, but argues against a role in cooperating alongside syt-1 in a pV_r-based facilitation mechanism.

Since dual fusion-sensor models and other models depending on changes in pV_r (see Discussion) are unlikely to be sufficient, we next investigated mechanisms involving an activity-dependent regulation of the number of participating release sites. For this we extended the single-sensor 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 inter-stimulus 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²⁺ would stabilize the RRP/release site occupation by slowing down unpriming (Figure 7A). This made the steady-state RRP size dependent on the resting Ca²⁺ concentration and the modest dependence of this on the extracellular Ca²⁺ resulted in RRP enlargement with increasing extracellular Ca²⁺ (Figure 7B), in agreement with recent findings on central synapses (Malagon et al., 2020). This model (like the dual fusion-sensor models depicted in Figure 6 and Figure 6—figure supplement 1) includes two different Ca²⁺ sensors, but the major difference is that these Ca²⁺ sensors operate to regulate two separate sequential steps (priming and fusion). Indeed, this scenario aligns with reports of a syt-7 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 paired-pulse experiment for low and high extracellular Ca²⁺ 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₁ priming is submaximal (~41%) for 0.75 mM extracellular Ca²⁺, but near complete (~99%) at 10 mM extracellular Ca²⁺. At low extracellular Ca²⁺ the elevation of Ca²⁺ caused by AP₁ results in a sizable inhibition of unpriming, leading to an increase (‘overfilling’) of the RRP during the inter-stimulus interval. With this, more primed SVs are available for AP₂, causing facilitation (green line in Figure 7C). In contrast, at high extracellular Ca²⁺ 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²⁺ concentration, AP₁ induces a larger Ca²⁺ 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²⁺ independent rate k_rep). Because Ca²⁺ acts in RRP stabilization, not in stimulating forward priming, this model (unlike the dual fusion-sensor models in Figure 6 and Figure 6—figure supplement 1) did not yield asynchronous release in-between 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²⁺.

Figure 7 with 3 supplements see all

Download asset Open asset

Figure 7—source data 1

Simulation data for graphs in Figure 7E, G, I; Figure 7—figure supplement 1B, D, E; Figure 7—figure supplement 2D, F, H and Figure 7—figure supplement 3.

https://cdn.elifesciences.org/articles/51032/elife-51032-fig7-data1-v2.xlsx

Download elife-51032-fig7-data1-v2.xlsx

In this model we assumed a Ca²⁺ 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 fusion-sensor models), together with K_M,prim, the Ca²⁺ affinity of the priming sensor, and u, its Ca²⁺ cooperativity. These values together define the Ca²⁺-dependent unpriming rate (see Table 2 for best fit parameters). The total number of fitted parameters (5) was the same as for the dual fusion-sensor models (Figure 6 and Figure 6—figure supplement 1). Figure 7D–I present the results. It is clear that both eEJC₁ amplitudes and PPR values were described very well with this model at all extracellular Ca²⁺ concentrations. In addition, the variance-mean relationship of the eEJC₁ was reproduced satisfactorily, except for a small variance overshoot for the highest extracellular Ca²⁺ concentrations (Figure 7I, see Discussion). Fitting of the unpriming model with a Ca²⁺ 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 time-dependence of the facilitation by simulating PPR values for various inter-stimulus intervals at different extracellular Ca²⁺ 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₁ 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 fusion-sensor model, PPR values were fairly robust to changes in Ca²⁺ influx (note the different scales on Figure 7J,K and Figure 6J,K). Moreover, because this mechanism does not affect the Ca²⁺ sensitivity of SV fusion, facilitation was achieved without inducing asynchronous release (Figure 7D).

We also investigated an alternative model based on Ca²⁺-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²⁺-dependent. In order to avoid site activation during AP₁, 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²⁺-independent reaction. This could mean that priming occurs in two-steps, with the first step being Ca²⁺-dependent. Similar to the unpriming model presented above, the modest increase of intracellular Ca²⁺ with extracellular Ca²⁺ 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²⁺-dependent eEJC₁ 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.

Why do n_site/priming-based 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²⁺ 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 paired-pulse experiment (0.75 mM extracellular Ca²⁺) in greater detail (Figure 8). Panel 8A, similarly to Figure 5A, shows example stochastic simulations (at external Ca²⁺ concentration 0.75 mM, to illustrate facilitation). The best fit parameters of the unpriming model predicted a larger Ca²⁺ influx (1.64-fold and 3.05-fold larger Q_max value) than the single- and dual fusion-sensor models (Table 2). The larger Ca²⁺ 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₁, vesicles prime into empty sites across the entire distribution, allowing AP₂ to draw again from the entire distribution. During this time, the increased residual Ca²⁺ causes overfilling of the RRP, that is more release sites are now occupied, giving rise to more release during AP₂. Notably, the AP₂-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 fusion-sensor model.

Figure 8

Download asset Open asset

Flies were kept under standard laboratory conditions as described previously (Sigrist et al., 2003) and reared on semi-defined medium (Bloomington recipe) at 25°C, except for GCaMP6m and synapGCaMP6f flies which were kept at room temperature, and Ok6-GAL4/+ (Figure 2, Figure 2—figure supplement 1, Figure 4 panel B-E and G, Figure 6 panel D-G 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 Ok6-GAL4/+ (Ok6-Gal4/II crossed to w[1118]; panel D-G, (I). Figure 7—figure supplement 1: Ok6-GAL4/+ (Ok6-Gal4/II crossed to w[1118]). Figure 7—figure supplement 3: Ok6-GAL4/+ (Ok6-Gal4/II crossed to w[1118]; panel C-F, (H).

The following stocks were used: Ok6-GAL4/II (Aberle et al., 2002), UAS-Unc13A-GFP/III (Böhme et al., 2016), elav-Gal4/I (Lin and Goodman, 1994). The following stock were obtained from the Bloomington Drosophila Stock Center: P{w[+mC]=Mhc-SynapGCaMP6f}3–5/III (Newman et al., 2017) and w[1118]; P{y[+t7.7] w[+mC]=20XUAS-IVS-GCaMP6m}attP40. The following stock was obtained from Kyoto Stock Center: P84200/IV.

Sample preparation, EM image acquisition and the quantification of docked SV distances to the AZ center (center of the electron dense ‘T-bar’) are described in Böhme et al. (2016); Reddy-Alla et al. (2017). The Rayleigh distributions were fit to the distances of docked SVs to the T-bar pedestal center, which had been collected in two EM datasets; analyses of these datasets were published in two previous studies, (Reddy-Alla 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.

The distances between Ca²⁺ 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 cross-section of the active zone, we integrate this distribution around a circle to obtain the two-dimensional 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:

In 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)$.

Sample 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. Third-instar w[1118] larvae were put on a dissection plate with both ends fixed by fine pins. Larvae were then covered by 50 µl of ice-cold hemolymph-like saline solution (HL3, pH adjusted to 7.2 [Stewart et al., 1994]: 70 mM NaCl, 5 mM KCl, 20 mM MgCl₂, 10 mM NaHCO₃, 5 mM Trehalose, 115 mM D-Saccharose, 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 ice-cold HL3 and then fixed for 5 min with ice-cold methanol. After fixation, samples were briefly rinsed with HL3 and then blocked for 1 hr in 5% native goat serum (NGS; Sigma-Aldrich, MO, USA, S2007) diluted in phosphate buffered saline (Carl Roth Germany) with 0.05% Triton-X100 (PBT). Subsequently dissected samples were incubated with primary antibodies (guinea-pig 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 fluorescence-labeled secondary antibodies (goat anti-guinea 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 (Max-Planck Institute for Biophysical Chemistry, Group of Stefan Hell) on high-precision glass coverslips (Roth, Germany, LH24.1). Two-color STED images were recorded on a custom-built STED-microscope (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⁻¹¹. The point spread function (PSF for deconvolution was generated using a 2D Lorentz function with its half-width and half-length fitted to the half-width and half-length 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; Wild-type 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 Unc13A-stained 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,y-coordinates 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 n-th 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 y-direction.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 y-axis of the histogram show percentage of the total amount of intensity.

Stage L3 larvae (n = 17; genotype: w[1118]; P{w[+mC]=Mhc-SynapGCaMP6f}3–5, Bloomington #67739) were fixed in ice-cold Methanol for 7 min and IHC-stained for BRP (mouse anti-Nc82, 1:1000; secondary AB: goat anti-mouse Cy5 1:500). Confocal images of the preparations were taken and processed as described in Reddy-Alla 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 ImageJ-function ‘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,y-coordinate to all others were determined using the MATLAB function pdist2, resulting in a square matrix containing all possible inter-AZ 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 1-NND through 4-NND in Figure 3A. The mean distance of the four nearest neighbouring AZs (1–4-NND) 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).

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 haemolymph-like solution without Ca²⁺ (HL3, pH adjusted to 7.2 Stewart et al., 1994: 70 mM NaCl, 5 mM KCl, 20 mM MgCl₂, 10 mM NaHCO₃, 5 mM Trehalose, 115 mM D-Saccharose, 5 mM HEPES) on Sylgard (184, Dow Corning, Midland, MI, USA) and transferred into the recording chamber containing 2 ml of HL3 with CaCl₂ (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 P-97 micropipette puller (Sutter Instrument, CA, USA) and filled with 3 mM KCl. Signals were low-pass 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.

Because the presynaptic terminals of the Drosophila larval NMJ are not readily accessible to electrical recordings of Ca²⁺ currents, the saturation behaviour of Ca²⁺ influx as a function of extracellular Ca²⁺ concentrations was measured. We did so by engaging the fluorescent Ca²⁺ indicator GCaMP6m (Genotype: w[1118]; P{y[+t7.7] w[+mC]=20XUAS-IVS-GCaMP6m}attP40, Flybase ID: FBti0151346), which we expressed presynaptically using OK6-Gal4 as a motoneuron-specific driver. Third instar larvae heterozygously expressing the indicator were used in experiments as follows. Dissection took place in Ca²⁺-free, standard hemolymph-like solution HL-3 (in mM: NaCl 70, KCl 5, MgCl₂ 20, NaHCO₃ 10, Trehalose 5, Sucrose 115, HEPES 5, pH adjusted to 7.2) (Stewart et al., 1994). After dissection on a Sylgard-184 (Dow-Corning) block, larvae were transferred to the recording chamber containing HL-3 at varying CaCl₂ concentrations (see below). The efferent motoneuron axons were sucked into a polished glass electrode containing a chlorided silver-wire, which could be controlled via a mechanical micromanipulator (Narishige NMN25) and was connected to a pipette holder (PPH-1P-BNC, NPI electronics) via a patch electrode holder (NPI electronics), and connected to an S48 stimulator (Grass Technologies). Larvae were then recorded using a white-light source (Sutter DG-4, 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://micro-manager.org) on an upright microscope (Olympus BX51WI) with a 60x water-immersion 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²⁺ concentration (0.75 mM) and then repeated in the same larva at increasing Ca²⁺ concentrations (in mM 1.5, 3, 6) by exchanging the extracellular solution. To achieve a situation with no Ca²⁺ influx, a final recording was conducted where the bath contained HL-3 without CaCl₂ and instead 8.3 mM EGTA (this solution was made by diluting 2.5 ml of a 50 mM stock solution in H₂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²⁺-containing solutions by adding 2.5 ml H₂O to 12.5 ml of HL3 before CaCl₂ was added at above mentioned concentrations.

Analysis of 5 Drosophila 3^rd instar Larvae was done after automated stabilization of x,y-movement in the recordings (8-bit multipage .TIF-stacks, converted from 16 bit) as described previously (Reddy-Alla 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 (frame-wise) from the signal, separately for each single recording. The quantification was then performed individually for each Ca²⁺ 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²⁺]_ext is the extracellular Ca²⁺ concentration. K_M,fluo (best fit value: 2.679 mM) is the concentration of extracellular Ca²⁺ at which fluorescence was half of F_end. The exponent m indicates a cooperative effect of the extracellular Ca²⁺ concentration on the fluorescence increase, which was constrained to a value of 2.43 (unitless) based on the described Ca²⁺ 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).

We 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₁ and A₂ represent the amplitudes of the first and second release, respectively, capital ‘E’ denotes the mean of a stochastic variable (e.g. EA₁), and a₁ and a₂ 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₁ = a₁ 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 non-linear function of two random variables is not equal to the non-linear function evaluated in the means). An example is shown in Figure 4—figure supplement 1, where the single-sensor model was simulated with varying amounts of Ca²⁺ influx (by varying Q_max). The most left blue point, for example, is significantly higher than the deterministic estimate (p=4e-16, one-sample t-test). This motivates the use of stochastic simulations for correct estimation of the PPR.

All 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.

1. Given a set of parameters, we first ran deterministic Ca²⁺ simulations in space and time in the presynapse at the desired extracellular Ca²⁺ concentrations.

2. 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²⁺ concentrations for the exocytosis simulation.

3. The simulation of the models for Ca²⁺ binding and exocytosis was performed for each SV position with the Ca²⁺ transients giving rise to the changing reaction rates.

4. 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.

Simulation of Ca²⁺ 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²⁺ 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²⁺ concentrations were simulated in space and time in a cylinder-shaped 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²⁺ from adjacent AZs (Meinrenken et al., 2002) and a volume-distributed 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². 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³. Increasing the height further had no effect on the Ca²⁺ transients.

The total amount of charge flowing into the cell was assumed to relate to extracellular Ca²⁺ in a Michaelis-Menten-like 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²⁺ 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 ${\displaystyle \sigma =\frac{0.360}{2\sqrt{2\cdot \mathrm{l}\mathrm{n}\left(2\right)}}=\mathrm{0.153.}}$

The CalC simulation output were data files that contained the spatio-temporal intracellular Ca²⁺ 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,y-plane and at time points with MATLAB’s built-in interpolate functions when computing the reaction rates of the system at a given time point.

The resting Ca²⁺ concentration was assumed to relate to the extracellular Ca²⁺ concentration in a similar way as during stimulation, such that

with ${\left[C{a}^{2+}\right]}_{max}=190nM$

For designation and value of Ca²⁺ parameters, see Table 1.

In 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²⁺ 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.

The models are summarized in Figures 4A, 6A and 7A, and Figure 7—figure supplement 3A,B. In the following equations the single-sensor, dual fusion-sensor, and unpriming models are all described. The site activation model is a combination of the equations for the single-sensor model and the site activation equations described below. The red text denotes terms that are unique to the dual fusion-sensor model, blue text 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 single-sensor model, dual fusion-sensor model and unpriming model: