A sequential two-step priming scheme reproduces diversity in synaptic strength and short-term plasticity

Fig. 2 A and B illustrate eEPSCs and the time courses of mean quantal content (m) in response to stimulus trains consisting of 15 APs (0.5, 1, and 2 Hz) or 40 APs (5–200 Hz). For a given synapse, m estimates for the initial eEPSCs (m₁) were similar across all f_stim. However, the mean m₁ across all f_stim varied nearly 10-fold between synapses (77–739 SVs, coefficient of variation (CV) = 0.42; Fig. 2D). The average over all 35 mean m₁ values was 377 ± 27 SV (Fig. 2B). During trains, the average m_j decreased monotonically toward a depressed steady state for all but the highest f_stim. During 200 Hz stimulation, net facilitation was observed for the average response—i.e., the PPR (m₂/m₁) was on average >1 (SI Appendix, Fig. S1A). PPR, when plotted against ISI, converts from facilitation into depression with a time constant of ∼21 ms (SI Appendix, Fig. S1A) (48). Plotting PPR_{200 Hz} as a function of m₁ for individual synapses revealed large diversity and a negative correlation—i.e., synapses with larger m₁ predominantly showed paired-pulse depression while those with smaller m₁ often exhibited PPF (Fig. 2D). Such correlation is usually interpreted as an indication of heterogeneous initial p_fusion (p_fusion,1), because strong depression in synapses with high p_fusion,1 would occlude PPF (49). The two-step priming scheme (Fig. 1 A and B) allows an alternative view, as detailed in SI Appendix: It provides an approximate expression for calculating p_fusion,1 from PPR and STD during 10 Hz trains (Eq. 31) for each synapse. Plotting such p_fusion,1 estimates versus the respective m₁ values revealed only weak correlation (Fig. 2E). The mean p_fusion,1 amounted to 0.43 with a CV of 0.20. The latter was substantially smaller than the CV of m₁ values, indicating that the variability of p_fusion,1 is not the principal cause for heterogeneous synaptic strength. It is rather due to a variable size of the subpool of fusion-competent SVs at rest (SP_TS,rest), which is readily calculated for each synapse as the ratio m₁/p_fusion,1. Plotting SP_TS,rest estimates as a function of the respective m₁ values revealed a strong linear correlation with an expected slope of ∼1/0.43 (Fig. 2F). The estimated mean SP_TS,rest (880 ± 57 SVs) had a CV of 0.38, which was close to the CV of m₁ values.

In sum, a large heterogeneity among synapses with respect to their SP_TS size at rest can explain strong variability in m₁, while p_fusion varies only little.

We next derived suitable initial guesses for model parameters from NTF analysis. NTF decomposes complex data sets into a linear combination of components (50) and was recently adapted for the analysis of eEPSC trains (51). The algorithm considers time courses of quantal contents during trains as superpositions of contributions by two or more types of signal sources and assumes that their differential relative contributions account for STP diversity among synapses. In the context of the model (Fig. 1A), these sources correspond to the release contributions by SVs residing in certain states (such as TS and LS) prior to stimulation. NTF analysis returns the time courses of individual contributions as BFs (51). For a given f_stim, BFs are the same for all synapses. They are normalized to a cumulative sum of 1, such that the product of a BF and the corresponding synapse-specific train quantal content (M) represents the time-resolved quantal release of that component for a given synapse (Fig. 3B and SI Appendix, Fig. S2 A and B). Two-component NTF fits provide BFs for release contributed by preexisting SV_TSs (BF_TS; Fig. 3A) and for the combined remaining release originating from SVs, which had been either loosely docked at stimulation onset or were newly recruited during the train (BF_LS,RS; Fig. 3B). The latter release component can be decomposed by subsequent 3-component NTF analysis (Fig. 3C and SI Appendix, Fig. S2 A and B). In short, NTF provides time courses of release contributions of SVs or sites, which had been in one of the states of the model (ES, TS, and LS; Fig. 1A) at stimulation onset.

NTF analysis cannot provide initial guesses for all model parameters, but it is very instrumental in constraining some of them. Two-component NTF, which is robust (51), separates release contributed by preexisting SV_TSs from other contributions. BF_TSs decay rapidly and approach zero after ∼5 APs (Fig. 3 A). Provided that all preexisting SV_TSs are consumed, their initial value represents p_fusion,1 (0.39 ± 0.004 when averaged over all six f_stim), and the train quantal content (M_TS) associated with BF_TS represents an estimate for the subpool SP_TS at rest (SP_TS,rest; average 961 ± 68 SVs). BF_TSs for f_stim of 5 to 20 Hz are strikingly similar when plotted against stimulus number (Fig. 3 A, i). They decay exponentially, indicating nearly constant p_fusion throughout trains at these frequencies. BF_TSs for ≥50 Hz deviate from this pattern (Fig. 3 A, ii). Their second values are larger, followed by a steeper decline, which is especially prominent for 200 Hz. Quantitative analysis of the time courses of p_fusion during trains as derived from BF_TS (51) indicates a small decrease in p_fusion for f_stim of 5 to 20 Hz, while for ≥50 Hz p_fusion increases, consistent with the observed net facilitation (Fig. 3D). Fig. 3B shows the time courses for quantal release contributed by those SVs, which were not in the TS state at stimulus onset (Release_LS,RS). Again, for 5 to 20 Hz, these time courses are strikingly similar when plotted against the stimulus number, indicating similar contributions to release regardless of ISI duration (Fig. 3 B, i).

Three-component NTF decomposition allows the separation of release contributions from preexisting SV_LSs (Fig. 3 C) and newly recruited SVs (SI Appendix, Fig. S2A) (51). Provided that all preexisting SV_LSs are consumed, the train quantal content (M_LS) associated with the BF_LS represents an estimate for SP_LS at rest (SP_LS,rest; average 1,078 ± 75 SVs). BF_LSs start at a value very close to zero since preexisting SV_LSs cannot fuse during the first AP but first need to undergo the LS → TS transition. For all frequencies, BF_LSs quickly increase to a maximum value after 2 to 3 APs before decaying exponentially to near zero. Again, the time courses of BF_LSs for 5 to 20 Hz are nearly indistinguishable (Fig. 3 C, i). As pointed out previously (51) and formally proven here (Eqs. 10–20), this finding is consistent with a transition of a constant fraction (s₂) of SP_LS to SP_TS subsequent to an AP. The average BF_LS for 5 to 20 Hz trains can be used to obtain an estimate for s₂ (0.11; Fig. 3 C, i, Inset). The BF_LSs for 50 to 200 Hz reach higher peak values and thereafter approach zero with faster decay kinetics (Fig. 3 C, ii).

At steady state, the loss from SP_TS due to SV fusion is compensated by an equal number of SVs replenishing SP_TS. Because the time course of BF_LSs (Fig. 3 C, i) along with that of the release contributed by newly recruited SVs (Release_RS; SI Appendix, Fig. S2 A, i) are frequency invariant for f_stim = 5 to 20 Hz, subpool occupancies approach very similar steady-state values regardless of ISI durations. That way, SVs are supplied to SP_TS at a rate linearly increasing with f_stim, which yields a frequency invariant steady-state quantal content (m_ss) at 5 to 20 Hz (SI Appendix, Fig. S1 B–D) (52, 53). We refer to this feature as balanced priming.

In summary, NTF analysis provides initial estimates for key model parameters, such as initial p_fusion, the size of SP_TS and SP_LS at rest, and the fraction of SP_LS transferred to SP_TS subsequent to each AP. Several more parameter estimates can be derived from analytical expressions regarding model predictions for, e.g., PPR and steady-state release at low-frequency stimulation (Eqs. 22–25). The time courses of the NTF components and their frequency invariance at 5 to 20 Hz suggests that APs act autonomously in this frequency range by releasing a certain fraction of SVs (p_fusion) and shifting the same amount between subpools, irrespective of the interval between successive stimuli.

After initializing model parameters (Fig. 1B) as described, numerical simulations reproduce steady-state depression (D_m = m_ss/m₁) for f_stim = 1 to 20 Hz quite well (Fig. 4D, dotted trace). We assumed for both priming rate constants k₁ and k₂ fixed resting values (k_1,rest and k_2,rest) and linear slopes σ₁ and σ₂ to describe their Ca²⁺ dependence (Eqs. 4, 5). The resting sizes of subpools, SP_LS,rest and SP_TS,rest, are thus determined by k_1,rest, b₁, k_2,rest, and b₂ (Eqs. 1–5). During stimulation, SP_TS partially depletes, such that release decreases toward a steady-state value m_ss for which SV consumption and replenishment are at equilibrium (Fig. 4D). At f_stim ≤2 Hz, this balance depends on all the parameters listed here, and increasing f_stim results in reduced replenishment and lower m_ss due to shorter ISIs (Fig. 4D). However, at f_stim of 5 to 20 Hz, the Ca²⁺-dependent increments of k₁ and k₂ above their resting values dominate the priming rates. Unpriming, determined by b₁ and b₂, becomes negligible, and the net movement of SVs along the kinetic scheme occurs nearly exclusively in the priming direction. The resulting steady-state occupancy of subpools is dominated by p_fusion and the integrals s₁ and s₂ over priming rates between consecutive APs (Eqs. 10, 12, and 24), the latter representing fractions of upstream subpools converted per AP to the corresponding downstream subpools. s₁ and s₂ are related to the slopes σ₁ and σ₂ according to Eq. 12. They are independent of f_stim, as long as [Ca²⁺] transients are short relative to ISIs and [Ca²⁺]-dependent priming dominates. Thus, m_ss tends toward a nearly frequency-independent plateau value (Eq. 24) since both SV consumption and replenishment increase approximately linearly with frequency (Fig. 4D, dotted trace) (52, 53). Recovery from STD following conditioning 10 Hz stimulation is well approximated by single exponentials with time constants ∼4 s (54), which can be reproduced in simulations by appropriate selection of k_1,rest and k_2,rest.

In sum, STP during 1 to 20 Hz eEPSC trains is well described by a balanced priming scheme, as suggested by NTF, in which each AP transfers constant fractions (s₁ and s₂) of SVs from the respective upstream subpools to SP_LS and SP_TS and releases an almost constant fraction of SP_TS (= p_fusion). This leads to increasing STD (decrease of m_ss) for increasing f_stim up to 5 Hz, when Ca²⁺-independent rates and Ca²⁺-dependent rates are of similar magnitude. However, for f_stim ≥5 Hz, when Ca²⁺-dependent rates dominate, occupancy of both SP_TS and SP_LS is nearly frequency-independent, leading to relatively stable m_ss despite increasing f_stim. For f_stim >20 Hz, however, this simple scheme will have to be extended, as detailed below.

For f_stim ≥50 Hz, large deviations are evident between experimental data and the basic model, which therefore needs to be extended to reproduce STP for f_stim up to 200 Hz. Inspection of the respective BFs reveals several changes at high f_stim: 1) p_fusion increases during trains, which is reflected in elevated second and third values of BF_TS, followed by a more rapid decay (Fig. 3 A, ii and D); 2) release contributed by those SVs not residing in TS prior to stimulation develops a peak around the third to fifth AP before decaying to steady-state levels lower than those at ≤20 Hz (Fig. 3 B, ii and C, ii); and 3) steady-state release contributed by newly recruited SVs decreases with increasing f_stim, causing stronger steady-state depression (Figs. 3 B, ii and 4D and SI Appendix, Fig. S2 A, ii). These changes in synapse behavior, which are likely caused by a summation of [Ca²⁺] transients or an incomplete relaxation of internal states of the release machinery, prompted us to consider options for extending the basic sequential model (Fig. 1). This also required slight adjustments of the parameters used so far due to some overlap of high- and low-frequency features.

The increase in p_fusion is likely a consequence of Ca²⁺ current facilitation (55–57), which we simulated with an empirical model (Eqs. 37–41) to reproduce NTF analysis (Fig. 3D). However, the increase in p_fusion alone is insufficient to account for the measured PPF (48). In addition, the early peak at 200 Hz of the release component, reflecting the contribution by SVs that were not tightly docked prior to stimulation (Release_LS,RS; Fig. 3 B, ii), is not reproduced by the basic sequential model. We therefore explored additional mechanisms potentially accounting for release facilitation and enhanced priming at higher f_stim.

For f_stim ≥50 Hz, AP-induced global [Ca²⁺] transients are expected to summate (Fig. 4C). Hence, a supralinear relationship between priming and [Ca²⁺] might generate release facilitation. However, simulations agreed only poorly with experimental data when we implemented this mechanism. Considering recent ultrastructural data of presynaptic active zones (AZs), which document rapidly reversible priming (37), we felt compelled to introduce a labile tightly docked SV state (TSL) and corresponding SV_TSLs. We assume that each AP converts a certain fraction κ of SP_LS into such labile SV_TSLs, which constitute SP_TSL (Eqs. 43–46). They have the same p_fusion as SV_TSs, but their backward transition LS ← TSL (b₃) is more rapid than that of LS ← TS (b₂). Thus, they contribute “extra” release to all responses of a high-frequency train except eEPSC₁, thereby enhancing PPR. They do not contribute to release, when ISIs are much longer than their mean decay time ($1/b_3$). A labile priming state can be implemented in series between LS and TS or as a parallel branch to the LS → TS transition. For both options, model parameters can be found, which reproduce STP accurately, even when the p_fusion of SV_TSs and that of SV_TSLs are constrained to have the same value. For brevity, we describe here only the version with the parallel branch (Fig. 1C). SP_TSL and SP_LS are incremented and decremented, respectively, following each AP by an amount of κ · SP_LS. During ISIs, the time courses of subpool occupancies including SP_TSL are described by the differential equations Eqs. 43–46. By introducing this extension of a TSL together with the facilitation of p_fusion, we were able to numerically simulate release during high-frequency trains with high fidelity (Fig. 5B) and identified κ = 16% and b₃ = 11.1 s⁻¹ as best parameters (SI Appendix, Table S2). These numbers are compatible with recent EM data showing that the number of SVs tightly docked at the AZ is transiently increased 5 ms after an AP but has relaxed back to normal when flash-freezing is initiated 100 ms later (37). Similar shifts between loosely and tightly docked SVs were reported to be associated with the induction of long-term potentiation (LTP) in hippocampal synapses (39) as well as with the beta-adrenergic modulation of parallel fiber LTP in the cerebellum (40).

In sum, release facilitation during the onset of high-frequency trains is generated by an increase in p_fusion, by a more rapidly replenished SP_TS due to accelerated priming when global [Ca²⁺] summates, and further by a transient filling-up of SP_TSL.

At f_stim ≥50 Hz, a progressive decline of release after the second or third stimulus toward a lower m_ss is observed (Figs. 4 A, ii and 5B). It results in a saturation of the relationship between release rate and f_stim (SI Appendix, Fig. S1C). We explored two options for modeling this: 1) a saturation of the priming rate with increasing f_stim and 2) a delayed availability of recently used sites for SV docking by introduction of a refractory release site state (58, 59).

As the simplest case, we describe here the replacement of the linear relationship between k₁ and [Ca²⁺] (Eq. 4) by a Michaelis–Menten (MM) type saturating relationship (Eq. 42). This introduces a parameter K_0.5, which is the [Ca²⁺] at half-maximum k₁. A K_0.5 of 280 nM and a slight adjustment of other parameters adequately reproduces experimental data for all f_stim (Fig. 4A and SI Appendix, Table S2). Note that the simulated mean steady-state depression D_m = m_ss/m₁ exactly superimposes on experimental data over two orders of magnitude of f_stim (Fig. 4D).

Having established an adequate fit to the experimental data, we then examined time courses of model quantities such as effective [Ca²⁺], p_fusion, and subpool occupancies (Figs. 4 B and C and 3D). At 10 Hz stimulation, the occupancy of SP_LS remains relatively constant during the train, while SP_TS rapidly depletes to a level of 0.310 times its resting value. This decay is the major cause of observed STD at 10 Hz (m_ss/m₁ = 0.312 ± 0.013). At 200 Hz stimulation, both SP_LS and SP_TS strongly deplete, and the depletion of the former is a consequence of the MM type saturation of k₁. Depletion of subpools is partially compensated by an increase in p_fusion (Fig. 3D). Model predictions for steady-state subpool sizes and depression are plotted against f_stim in Fig. 4 D and E, illustrating that the occupancy of SP_LS at steady state remains nearly constant and close to its initial value (SP_LS,rest/N_tot = 0.44; SI Appendix, Table S2) for all f_stim ≤10 Hz (Fig. 4D). At 200 Hz stimulation, ∼92% of the release sites are vacant at steady state (Fig. 4E). Simulated BFs are shown in SI Appendix, Fig. S3. Considering that NTF-derived BFs only provided initial guesses for various model parameters and that the latter were subsequently adjusted to optimize model predictions, the agreement between simulated (SI Appendix, Fig. S3) and NTF-derived (Fig. 3) BFs is reassuring.

Contributions of simulated NTF components to total release at 200 Hz are shown in Fig. 5D for two cases of differing initial subpool sizes. In comparison to the NTF decomposition of the measured mean 200 Hz eEPSC train (SI Appendix, Fig. S2 B, ii), Fig. 5D illustrates the differential contribution to release by SV_LSs following conversion into SV_TSs or SV_TSLs. In the two modeled synapses, this generates either pronounced net facilitation (Fig. 5 D, i) or faster and deeper depression (Fig. 5 D, ii) during 200 Hz trains.

As an alternative to a saturation of k₁, we introduced an empty but refractory release site state (ERS) (59–61). In this case, release sites vacated following SV fusion are not instantaneously converted into empty sites available for SV docking but shift into ERS, from which they become available for SV docking with a rate constant b₄ (Fig. 1 C, ii). A b₄ of 3.6 s⁻¹ and a slight optimization of other parameters (SI Appendix, Table S2) resulted in model predictions, which for our standard parameter set were hardly discernible from those of the alternative approach of using an MM–type saturation of k₁ (Fig. 1 C, i). Significant differences between the two types of models were only observed when the occupancy of release sites at rest was reduced below 70%.

In sum, decreasing steady-state release at f_stim ≥50 Hz indicates a saturation of the pool replenishment process, which can be simulated equally well by a saturation of the Ca²⁺-dependent acceleration of k₁ or else by introducing a refractory state of release sites, which after their use become available for new SV docking with first-order kinetics.

We next studied STP during 100 and 200 Hz eEPSC trains preconditioned by two or four APs delivered at 10 Hz, which causes pronounced release facilitation (62). Two examples for eEPSCs elicited by such a stimulus pattern, which are part of our standard NTF protocol, are shown in Fig. 6 A and B. Nearly all synapses showed more or less pronounced facilitation following 10 Hz preconditioning (Fig. 6C). The conditioning low-frequency stimulation depletes SP_TS, thereby exposing release facilitation at the onset of the high-frequency train due to the increase in p_fusion and because of the rapid conversion of SV_LSs to fusion-competent SV_TSs and SV_TSLs. The agreement between model predictions and average data was remarkable (Fig. 6D), considering that model parameters were largely determined on the basis of regular trains without preconditioning.

Finally, we tested the adequacy of the model against general features of eEPSC recovery from STD. As described above, the slow recovery after 10 Hz trains (54) can be reproduced by appropriate selection of k_1,rest and k_2,rest. Following high-frequency stimulation, the recovery time course gains a rapid component (10, 63). After cessation of stimulation, three processes overlap: a fast drop of p_fusion from a facilitated value back to p_fusion,1, a rapid decrease of SP_TSL, and an increase in SP_TS due to accelerated priming while [Ca²⁺] is elevated. The net result may either be an accelerated recovery or a transient decrease of synaptic strength, depending on the relative magnitudes of these processes and their kinetics (63). Ca²⁺-accelerated priming is sensitive to the [Ca²⁺] buildup during trains, which is determined by the decay of individual [Ca²⁺] transients. At ≤20 Hz, when individual [Ca²⁺] transients do not overlap, the Ca²⁺-dependence of k₁ and k₂ is relevant only in terms of the integral over transients—i.e., large transients with short duration are as effective in promoting priming as small transients with correspondingly longer duration. At ≥50 Hz, however, [Ca²⁺] transients summate. During stimulus trains, the increased rate of priming is balanced by release. After stimulation, refilling of partially depleted subpools continues at an increased rate until [Ca²⁺] has decayed back to [Ca²⁺]_rest. Lengthening the decay of the [Ca²⁺] transient at the expense of its amplitude shifts a larger proportion of its priming-promoting effect into the early recovery time course. In simulations, the amount of accelerated eEPSC recovery can be adjusted to match experimental data by varying the value of the model parameter τ_Ca (SI Appendix, Table S2). An example for a simulation of a stimulation pattern probing recovery is given in SI Appendix, Fig. S4. It should be noted, though, that experimentally observed decay time courses of global [Ca²⁺] are often biphasic (63–65). More detailed simulation of the Ca²⁺ signal relevant for controlling Ca²⁺-dependent priming (i.e., effective [Ca²⁺]) may improve the accuracy of model predictions.

NTF-based decomposition of eEPSC trains into release components rests on the following four assumptions: 1) Contributions to release are nonnegative. 2) Release can be decomposed into distinct components, each of which represents the contribution by SVs that had been in a certain state prior to stimulation (LS, TS, or undocked). 3) Forward transition rates of SV priming and fusion strongly dominate over backward transition rates in the f_stim range used for NTF analysis (5–200 Hz). 4) Heterogeneity among synapses with respect to initial strength and STP characteristics is primarily due to differences in the relative abundance of SV_LSs and SV_TSs at rest—i.e., individual synapses are endowed with variable fractions of docked SVs equipped with a mature release machinery.

NTF decomposes release into components constrained only by nonnegativity. Therefore, a given NTF fit result should not be regarded as a unique solution but rather as one of many options consistent with both the mean eEPSC trains of all synapses examined and their variability. Provided that the above-mentioned assumptions hold, NTF analysis delivers good initial guesses for model parameters, which together with trial-and-error parameter optimization leads to a number of important conclusions discussed in the following sections.

When analyzing 5 to 20 Hz eEPSC trains, 1) we postulated a linear relationship between [Ca²⁺] and the priming rate constants k₁ and k₂ and 2) we assumed that the integral over the AP-induced [Ca²⁺] transient is constant over that f_stim range. This results in balanced priming, with the average SV recruitment rate and the average release rate both being proportional to f_stim at steady-state conditions. However, the same conclusions can be reached with much-less-specific assumptions. Any process by which an AP causes the release of a certain fraction of fusion-competent SVs and also the transition of a certain fraction of immature SVs to the fusion-competent state will reach a frequency-independent steady state. This would hold even if recruitment were nonlinearly dependent on [Ca²⁺] as long as [Ca²⁺] transients summated only little, such that each AP exerted its Ca²⁺-dependent effect independently. This is largely the case for f_stim ≤20 Hz. The choice of the exact time course of [Ca²⁺] transients has little influence on the outcome of simulations in this f_stim range, because different [Ca²⁺] transients will similarly influence release as long as they have the same time integral and decay within the ISI back to [Ca²⁺]_rest (Eqs. 12–14).

The simplicity of conditions in the 5 to 20 Hz range (almost constant p_fusion, balanced loss and gain for SP_LS, fixed fraction of source pools transferred to downstream pools per AP) allowed us to derive a simple equation for calculating p_fusion,1 from the two easily measurable quantities, PPR = m₂/m₁ and D_m = m_ss/m₁. In this expression (p_fusion,1 = (1 – PPR)/(1 – D_m); Eq. 31), the numerator recapitulates the standard use of PPR as an indicator for p_fusion (14, 66). The equation shows that 1 – PPR is actually a good estimate for p_fusion in case of complete steady-state depression (D_m = 0), provided that all simplification made in the derivation of Eq. 31 can be applied. However, SV pool depletion is usually far from complete for 5 to 20 Hz stimulation for which Eq. 31 holds. Thus, the correction by the denominator is substantial. The p_fusion,1 value obtained that way (0.43) is quite close to that derived from NTF analysis (0.39), the latter being used in all numerical simulations. Estimates for p_fusion,1 obtained according to Eq. 31 for individual synapses show less variability than corresponding m₁ values, which argues against p_fusion,1 being the main source of heterogeneity in synaptic strength among synapses (Fig. 2 D–F).

For a parallel priming and fusion scheme, the mean p_fusion is simply the weighted average of the p_fusion values pertaining to individual subpools of fusion-competent SVs. The mean p_fusion may be low, if “reluctantly releasing” SVs (9, 18, 67, 68) contribute to release. In the context of the sequential priming scheme discussed here, our analysis suggests a relatively high p_fusion for SV_TSs, which are the only fusion-competent SVs in resting synapses (Fig. 1B). In contrast, “traditional” estimates often calculate p_fusion as the ratio of the quantal content of a single synaptic response over the size of the pool of readily releasable SVs as determined by high-frequency stimulus trains (69). Under the conditions described here, such stimulus trains deplete not only SP_TS but rather the sum SP_TS + SP_LS due to the rapid LS → TS transition during high-frequency stimulation. Release probability as measured by traditional methods, therefore, is the product p_fusion · SP_TS/(SP_LS + SP_TS), i.e., p_fusion times the probability of a docked SV being in the tightly docked state equipped with a mature release machinery (51). Considering a sequential priming scheme with only one fusion-competent state, p_fusion as determined by NTF or model fitting is a quantity strictly reflecting the SV fusion process, whereas p_fusion as determined by traditional methods is a quantity depending on the fusion process and on the priming equilibria at rest (70, 71) (SI Appendix, Fig. S5). This insight, if applicable to a given type of synapse, has strong implications for the interpretation of perturbations of synapse function by mutagenesis or pharmacological tools: An observed change in p_fusion is generally interpreted as a change in the fusion machinery, involving SNARE proteins, synaptotagmins, and complexins, or else reflecting changes in the microdomain [Ca²⁺] signal, which depends on Ca²⁺ currents, Ca²⁺ buffers, and coupling distances between voltage-gated Ca²⁺ channels (VGCCs) and Ca²⁺ sensors for SV fusion. According to our interpretation, an observed change in a traditional p_fusion estimate may well reflect changes in the LS ↔ TS equilibrium at rest, possibly involving presynaptic proteins such as Munc13, DOC2, CAPS, and synaptotagmin7 (27, 30). A modulatory influence of second messengers ([Ca²⁺], diacylglycerol (DAG), phosphatidylinositol 4,5-bisphosphate (PIP₂)) may show up in traditional quantal analysis as a change in p_fusion while NTF analysis reports it as a shift in the state of priming at rest. More precisely, the ratio of traditional p_fusion estimates over NTF-derived p_fusion estimates equals SP_TS/(SP_LS + SP_TS). Provided that conditions can be found for a given type of synapse for which Eq. 31 applies, p_fusion,1 can be estimated from PPR and D_m measured at a suitable frequency, and SP_TS,rest is readily calculated as the ratio m_1/p_fusion,1. Assuming further that pool estimates obtained by analyzing depleting high-frequency eEPSC trains represent the sum SP_LS,rest + SP_TS,rest, then both SP_TS,rest and SP_LS,rest can be determined approximately from a small number of measurements without kinetic modeling or NTF analysis.

For the f_stim range of 5 to 20 Hz, NTF provides quite stringent constraints on model parameters and eEPSC trains can be modeled as independently occurring release events largely unaffected by the ISI length. In contrast, our findings regarding STP mechanisms at f_stim ≥50 Hz can be interpreted in alternative ways. Consecutive release events during 50 to 200 Hz stimulation may interact in various ways: through the summation of global [Ca²⁺] transients, through kinetic limitations such as the saturation of priming rates, and via synergies, e.g., a supralinear dependence of p_fusion on increasing local [Ca²⁺] due to VGCC facilitation. Together, these interactions shape synaptic facilitation and depression, which can be prominently observed during 200 Hz stimulation during which initial and transient net facilitation of eEPSCs is followed by strong depression toward a small m_ss (Fig. 4 A, ii). To illustrate how a given interaction influences the release time course, we compare in Fig. 5 A and B experimental data with model predictions. The basic model, which is sufficient for describing low frequency–induced STP, fails to reproduce the 200 Hz data (Fig. 5A, solid trace). Allowing p_fusion to vary as predicted by NTF analysis and reproduced by the empirical y–z formalism (Fig. 3D and Eq. 37) is insufficient to explain the extent of PPF (Fig. 5A, dotted trace): Not only is m₂ far below the experimental value, but m₃ and m₄ are even more so. Adding the contribution of the TSL state to release leads to an adequate fit up to m₈ (Fig. 5B, dotted trace). However, it predicts subsequent release to rebound due to enhanced priming by the gradual buildup of [Ca²⁺] during 200 Hz trains. Only the additional implementation of k₁ saturation reproduces the experimentally observed release time course over the entire stimulation period (Fig. 5B, solid trace). Alternatively, the experimental data can be simulated satisfactorily by the introduction of an ERS into the model (Fig. 1 C, ii and SI Appendix, Table S2).

Previous studies of calyx synapses using long step depolarizations or presynaptic Ca²⁺ uncaging for triggering SV fusion identified two kinetically distinct release components mediated by two SV pools differentially coupled to presynaptic VGCCs (72, 73): a slowly releasing pool (SRP) and a fast-releasing pool (FRP) (68, 74, 75). The FRP accounts for the majority of SVs fusing during AP-evoked release, while the contribution of the SRP is only minor (76). The FRP may thus be regarded as the calyx–synapse equivalent of the readily releasable pool as estimated with AP trains in other synapses. Notably, SP_LS and SP_TS are not congruent with SRP and FRP but rather represent a subdivision of the FRP.

In a sequential priming scheme as used here, kinetic components may not be readily linked to specific state transitions. In contrast, it is intuitively easy to consider two kinetic components as independent contributions of two SV populations. Their intuitive tangibility explains why parallel kinetic schemes are particularly attractive. Previously, we thus described some of the STP features, based on a smaller data set, in terms of a parallel model consisting of a rapidly releasing SV subpool (called superprimed SVs) and a slowly releasing one (called normally primed SVs) (16) in analogy to studies at cultured hippocampal synapses (77). Superprimed SVs share many properties with SV_TSs. Normally primed SVs in the context of a parallel kinetic scheme represent fusion-competent SVs, albeit with a low fusogenicity (16). However, in the framework of the sequential scheme, proposed here, only SV_TSs are considered fusion competent. Therefore, and in view of the morphological evidence for “loose” and “tight” coupling, we consciously avoided the previous terminology.

Release mediated by SV_TSs has many features in common with the release contributed by so-called preprimed SVs, postulated for glutamatergic synapses in the cortex (78). Likewise, the sequential scheme proposed for cerebellar parallel fiber-molecular–layer interneuron synapses (22, 23) is related to the sequential scheme favored here: Both assume reversible priming steps in sequence, postulate certain fractions of upstream SV pools being transferred to a downstream pool upon AP arrival, and separate contributions to release in terms of the state of SVs prior to stimulation.

The identity of the [Ca²⁺] signal that regulates priming (effective [Ca²⁺]; SI Appendix, Table S2) is not unequivocally established. For f_stim of 5 to 20 Hz, each AP shifts certain constant fractions of SVs from one pool to the next downstream pool, which suggests a local, short-lived action of the effective [Ca²⁺] transients. If Ca²⁺-dependent priming depends on Munc13, then a local [Ca²⁺] transient similar to that triggering SV fusion may be most relevant since Munc13 is an integral part of the AZ (79–81). On the other hand, the priming rate may saturate far below the peak [Ca²⁺] within such local domains if the respective Ca²⁺ sensors have high affinity. A complete description of Ca²⁺-dependent priming would have to consider endogenous Ca²⁺ buffers, which shape the time course of both the local domain [Ca²⁺] (82, 83) and the global volume–averaged [Ca²⁺] (64, 65, 84, 85). Here we chose to explore the influence of the [Ca²⁺] time course with the simplest possible model, characterized by an AP-induced increment in [Ca²⁺] proportional to the Ca²⁺ influx, followed by an exponential decay toward [Ca²⁺]_rest (Fig. 4C and SI Appendix, Table S2; Eq. 6).

Fig. 5 C and D illustrates results of two simulations during which all model parameters were identical except for b₂, which was either increased or decreased to obtain a relative size of SP_TS,rest of 20% (Fig. 5 C, i and D, i) or 74% (Fig. 5 C, ii and D, ii), respectively. We then systematically varied b₂ or k_2,rest (Fig. 5E) and plotted the ratios m₂/m₁ and m₃/m₁ for simulated 200 Hz eEPSC trains as a function of the relative size of SP_TS,rest (Fig. 5 E, i). We observed negative correlations, reminiscent of the relationship between PPR and the measured m₁ (Fig. 2D). A similar negative correlation was found between the simulated steady-state depression m_ss/m₁ and simulated SP_TS,rest (Fig. 5 E, ii), illustrating that the magnitude of depression strongly depends on the relative abundance of SV_TSs at rest. Fig. 5 E, ii also plots the absolute m_ss as a function of the relative size of SP_TS,rest for simulated 200 Hz trains. As expected, m_ss is independent of the distribution of subpool sizes at rest since the Ca²⁺-dependent parameters of the priming scheme dominate at high frequencies over the small rate constants, which set the distribution of states at rest.

Recent studies emphasize variations in coupling distances between presynaptic VGCCs and docked SVs as a predominant mechanism generating differences in p_fusion (for review, see reference 86). Our implementation of NTF analysis, on the other hand, postulates differences in the maturation state of primed SVs as the main source of intersynapse variability. It is well conceivable that a different NTF implementation, allowing several states or types of release sites to contribute to release, would come up with BFs for such states with different p_fusion, as one might expect for different coupling distances. However, our assignment of variability to differences in the maturation state of primed SVs is self-consistent in the sense that SI Appendix, Eq. 31, when applied to individual synapses, confirms that variable p_fusion contributes little to their heterogeneity in synaptic strength (Fig. 2E). Of course, this does not preclude that spatial coupling between Ca²⁺ entry and the Ca²⁺ sensor for SV fusion may be heterogeneous.

The identification of molecular mechanisms causing functional synaptic diversity has recently attracted great attention due to the recognition of its importance for maximizing the capacity of information processing of neuronal networks (87) and of findings regarding the modulation of heterogeneity by astrocytes (88). Activity-induced changes in synaptic weight such as the redistribution of synaptic efficacy as a consequence of the induction of LTP (89) have emerged as important elements in the description of the dynamics of neuronal networks (90, 91). However, the understanding of molecular mechanisms underlying such phenomena has lagged behind. Such redistribution is a basic property of our model, controlled by the distribution of SVs between primed states at rest and not involving a change in p_fusion. In Fig. 5C, we compare simulated time courses of the release of two synapses, one with a low ratio of SP_TS/SP_LS (Fig. 5 C, i) and the other with a high ratio of SP_TS/SP_LS (Fig. 5 C, ii) at rest. Despite the 4.5-fold difference in initial synaptic strength (m₁ = 150 vs. 684 SVs), the cumulative release during the first 20 stimuli at 200 Hz is quite similar (2,198 vs. 2,695 SVs; SI Appendix, Fig. S5B). This suggests changes in the LS ↔ TS equilibrium at rest and priming proteins, such as Munc13s, as the molecular basis of a redistribution of synaptic efficacy.

Juvenile, posthearing onset (P14–16) Wistar rats of either sex were used. All experiments complied with the German Protection of Animals Act and with the guidelines for the welfare of experimental animals issued by the European Communities Council Directive. Acute brainstem slices were prepared similarly as previously described (16). See SI Appendix for details.

Whole-cell patch-clamp recordings were made from principal neurons of the medial nucleus of trapezoid body (MNTB) at room temperature as previously described (16). See SI Appendix for details.

NTF of eEPSC trains was performed similarly as previously described (51). See SI Appendix for details.

We used a 2-step priming and fusion scheme (Fig. 1) to numerically simulate SV fusion and short-term plasticity. Details of the model and approximations for steady-state conditions that were used to constrain model parameters are given in SI Appendix.