The spontaneous quantal release in cat presynaptic neuromuscular terminals, reported by Boyd and Martin (1956a), was fairly reproduced by our model fed with an α = 0.62 s^–1 value, obtained as the inverse of the experimental 1.61 s time constant (τ) of the time interval distribution of miniature end plate potentials (mepp_s). An unexpectedly large β = 100α (λ = β/α = 100 coefficient) and a ρ = 1.0 s^–1 recycling rate constant contributed to produce 148 ± 2 mepp_s at a 1.40 ± 0.10 s^–1 frequency (n = 250 simulations), quite similar to the 143 mepp_s recorded at a 1.43 ± 0.88 s^–1 frequency in the original study (Figure 2A). The experimental distribution of the intervals between mepp_s was fitted by the function n = n_T(△t/τ)e^−t/_τ (Fatt and Katz, 1952), where n_T is the number of quanta released and △t = 0.5s is the bin size.

The mepp_s frequency (Figure 2B) was proportional to α and inversely proportional to β. Our best explanation to this result was that the large β value keeps a reduced pool of primed vesicles, therefore, reducing the probability of spontaneous fusion. Figure 2C compares simulations of cat and frog spontaneous release. The value of β = 0.62 s^–1 (λ = 100) that reproduced the 1.43 s^–1 mepp frequency in cat recordings at 37°C (Boyd and Martin, 1956a) quadruples the β = 15 s^–1 (λ = 50) coefficient that reproduced the 2.5 s^–1 mepp frequency commonly recorded from frog synapses at 20°C (see Fatt and Katz (1952)). The larger rate constant values in mammalian neuromuscular synapses may reflect the characteristic higher physiological temperature of mammalian tissues.

The previous results may be explained in the following way. First, spontaneous release reflects spontaneous P → F transitions, and second, the small probability of spontaneous release depends on the large λ coefficient, which maintains a small pool of primed vesicles at rest. Since a majority of experimental evidence used here proceeds from experiments in frog, the simulations that follow used the α = 0.3 and β = 50 values, unless otherwise indicated.

Simulations of frog experiments made under low probability conditions (del Castillo and Katz, 1954a; Katz and Miledi, 1968) allowed a further analysis on the contribution of β to quantal release. The λ coefficients of the D ⇌ pP (λ₁) and pP ⇌ P (λ₂) transitions were varied independently, whereas the α = 0.3 s^–1, ρ = 1.0, and τ_e = 0.15 ms remained fixed. The λ₁ = λ₂ = 50 values reproduced transmission, as seen in the central chart of Figure 5.

The value of λ₁ markedly influenced the number of quanta discharged per impulse. A small λ₁ = 5 (β = 5α; top panels in Figure 5) that decelerates the D ↼ pP transition extended the range of classes in the distribution, which deviated from the predictions of the Poisson equation (p ≤ 0.05). Even the largest λ₂ = 500 value tested failed to compensate for the effect of a reduced λ₁. By contrast, a large λ₁ = 500 value constrained the amplitude mepp distribution to a small-class range that was predicted by the Poisson distribution, regardless of λ₂ (bottom plots in Figure 5). However, it will be shown below that this result only applies to release on single impulses as the large λ₁ = 500 values failed to reproduce short-term plasticity. In spite of that, the results in this section underscore the essential contribution of the backward D ↼ pP transition to maintain a small resting pool of primed vesicles.

The way by which the sequence of kinetic transitions affects the balance from facilitation and depression in frog was analyzed with the alternative protocol by Betz (1970). Experiments with high extracellular calcium concentration enhanced the release probability, while curare blocked acetylcholine receptors to render subthreshold transmission. Test impulses with different lags unveiled the time-dependence of depression. The long τ_e = 1.5 ms and our other frog parameters reproduced again the experimental results. As shown in Figure 7A, a briefer τ_e = 0.5 ms increased facilitation and abolished depression. By contrast, a longer τ_e = 2.0 ms eliminated facilitation but kept depression.

The short-term plastic dynamics of transmission upon a conditioning train followed by a test pulse are plotted in Figure 10, following the experiment by Mallart and Martin (1968, see Figure 5). The fraction of vesicles in each state was normalized to N₀ = 10,000. At rest, ∼98% vesicles are docked and the remaining 2% are decreasingly distributed in preprimed and primed states. About 300 vesicles (3%) fuse on the first impulse, as estimated by Katz and Miledi (1979), at 6°C. Therefore, ∼66% of vesicles that fuse were primed, the remaining arriving from immature states. Arrival of a second impulse encounters an increased population of preprimed and primed vesicles, thus evoking facilitation plus additional forward transitions in immature vesicles. After the third conditioning pulse, ∼25% of the total vesicle pool has fused. Such large release along with the slow recycling (F/N₀ panel in Figure 10) depress the response to the test impulse (Figure 10).

The similar rate constants of the forward transitions in the fusion complex (Li et al., 2007; Chapman, 2008; Sudhof and Rothman, 2009) requires three additions to reproduce the whole dynamics of transmission. First, the calcium-dependence of every forward transition. Second, the calcium-independent backward D ↼ pP ↼ P transitions become synchronized by spontaneously following the highly synchronic forward transitions. Third, a minimum of four-transitions is necessary and sufficient to reproduce the whole dynamics of neuromuscular transmission studied here.

It is interesting to note that the pP state buffers the effects of having logarithmic differences between the numbers of D and P vesicles. In absence of such buffering, a three-state sequence such as that suggested for the calix of Held synapse (Neher and Brose, 2018) results in exceedingly large amounts of release per impulse (Figure 2D). However, with adequate numbers of vesicles and release probabilities, the three-state sequence may reproduce the characteristic depression in the calix of Held (for review see von Gersdorff and Borst, 2002; Neher and Brose, 2018).

The balance between the forward and backward transitions explains the frequency-dependent non-linear fluctuations of the quantal output during facilitation and depression. The sequential transitions in the fusion complex on an impulse increase largely the pool of primed vesicles after synchronous exocytosis, producing facilitation upon rapid arrival of another impulse. Vesicle priming after the impulse is predicted by the model from the decreased or increased fusion latencies when kinetic steps are reduced or increased, respectively (Figure 3D). A corollary to this observation is that the whole essence for facilitation is that the forward D ⇀ pP ⇀ P reactions continue after the synchronous release, producing a transient accumulation of newly primed vesicles. Without such possibility, transmission would be dominated by depression.

Our data suggest that one calcium sensor may produce fusion in all forms of release. This result seems to contradict the generally accepted contribution of at least two calcium sensors with different calcium affinities (Yamada and Zucker, 1992; Kamiya and Zucker, 1994; for review see Zucker and Regehr (2002); Sun et al. (2007)) and different forms of synaptotagmin controlling vesicle fusion in the neuromuscular junctions (Pang et al., 2006; Wen et al., 2010) and central synapses (Sudhof, 2013; Kaeser and Regehr, 2014; Kavalali, 2015; Volynski and Krishnakumar, 2018). However, our model predicts that the calcium sensors promoting each transition on an impulse contribute to modulate the dynamics of release.

Central synaptic vesicles seem to carry different types of synaptotagmin (Jahn and Südhof, 1994; Takamori et al., 2006). While fast synchronous release is produced by the activation of synaptotagmins 1, 2, or 9 (Chapman, 2002; Pang et al., 2006; Xu et al., 2007, for review see Kaeser and Regehr (2014), Neher and Brose (2018)), asynchronous release is supposed to depend predominantly on the high calcium affinity synaptotagmin 7 (Wen et al., 2010; Bacaj et al., 2013, 2015; Turecek and Regehr, 2018). Accordingly, theoretical models of transmission with two or three calcium sensors reproduce well the electrophysiological data (Goda and Stevens, 1994; Dutta Roy et al., 2014). For convenience, it is useful to focus this section by analyzing first the evidence concerning asynchronous release.

Evidence has long suggested that facilitation and asynchronous neuromuscular release rely on the exact same mechanism (Rahamimoff and Yaari, 1973; Zucker, 1996). Our simulations are consistent with this idea. The generation of asynchronous mepp_s using a reduced τ_e value to eliminate the residual calcium effect on release suggests that asynchronous release is an exacerbated version of spontaneous release with increased numbers of primed vesicles after a conditioning impulse. Other line of evidence suggests that synaptotagmin 7 drives asynchronous release (Wen et al., 2010; Bacaj et al., 2013; Turecek and Regehr, 2018), although evidence has also shown that the same vesicles may participate on both modes of release (Grigoryev and Zefirov, 2015). However, in neuromuscular junction of zebra fish, elimination of synaptotagmin 7 reduces but does not abolish asynchronous release (Wen et al., 2010). Therefore, both, spontaneous fusion and synaptotagmin 7-driven fusion may contribute to asynchronous release in the neuromuscular junction. The question is when does synaptotagmin 7 produce its effects. According to our simulations, synaptotagmin 7 may have its effects on the calcium-dependent maturation steps rather than producing vesicle fusion. Such statement is supported by diverse effects of synaptotagmin stabilizing the D state and to the maturation of the vesicle fusion complex (Reist et al., 1998; Loewen et al., 2006; Mohrmann et al., 2013; for review see Bowers and Reist (2020)).

The residual calcium hypothesis for paired pulse facilitation by Katz and Miledi (1968) and the third or fourth order calcium-dependence of release (Dodge and Rahamimoff, 1967; Smith et al., 1985; Augustine and Charlton, 1986) predict that low residual calcium levels activate high-affinity calcium sensors to produce supralinear vesicle fusion in facilitation (Zucker and Lara-Estrella, 1983; Yamada and Zucker, 1992; Van der Kloot and Molgó, 1993; Vyshedskiy and Lin, 1997; Zucker and Regehr, 2002; Ma et al., 2015). Our model suggests the possibility that the calcium sensors producing facilitation are those activating the D ⇀ pP ⇀ P transitions, which increase the pool of vesicles ready for release. Synaptotagmin 7 has emerged again as a candidate in central synapses (Sugita et al., 2001; Bacaj et al., 2013, 2015; Jackman and Regehr, 2017; Turecek and Regehr, 2018). However, as mentioned above synaptotagmin 7 may be acting on the early molecular transitions. Therefore, according to our model, fusion is produced by one calcium sensor, while the modulation of the number of vesicles that fuse depends on the action of the calcium sensors on the early transition states with synaptotaagmin 7 being one such sensors.

Electron tomography shows that from the moment of docking, the fusion complex has formed intimate boundaries with calcium channels (Harlow et al., 2001; Nagwaney et al., 2009; Szule et al., 2012). The interactions between fusion complex proteins and calcium channels have been analyzed in detail (for review see Catterall et al., 2013; Gandini and Zamponi, 2021). Such configuration may permit calcium sensors to catalyze every kinetic transition, as opposed to central synapses in which calcium channels may be separated from fusion complexes in immature vesicles (Neher, 2015). Other proteins thought to be involved in docking and priming such as RIM, Munc13, rabphilin, and Bassoon/Piccolo, have calcium-binding domains which may contribute to these transitions (Friedrich et al., 2010; Nishimune, 2012; Gundelfinger et al., 2016; Lai et al., 2017).

Our results confirm the essential role of vesicle recycling on depression and predict that backward transitions contribute to the amplitude of depression. Two or more recycling modes in the neuromuscular junction (Rizzoli and Betz, 2005) and central synapses (Wu and Borst, 1999; Sakaba and Neher, 2001; Schneggenburger et al., 2002) suggest equal numbers of recycling vesicle pools (for review see Alabi and Tsien, 2012). However, with a single recycling rate constant, our model reproduced convincingly the balance between facilitation and depression as studied by Betz (1970). However, we cannot exclude that the slow time constant of recycling in our model is masking faster events including some displaying a calcium-dependence (Sakaba and Neher, 2001).

The four-state kinetic model with six kinetic transitions shown in Figure 1 is the basis to analyze the collective behavior of a pool of 10,000 identical vesicles (Rizzoli and Betz, 2005). Six R_j transitions correspond to those in Figure 1, with j being a stochastic discrete variable with values j = 1, 2,…6, that correspond to each kinetic transition. Each transition occurs with an equal probability a_j(x). The term a_j(x)dt is the probability that an R_j transition will occur in an infinitesimal time interval t + dt, when the system is in a state X(t) = (D(t), pP(t), P(t), F(t)) = x. Each R_j transition is characterized by two quantities: One is the system state x = D(t),pP(t),P(t),F(t), which reflects the number of vesicles at each kinetic state. The second quantity is the vector V_j(v_Dj, v_pPj, v_Pj, v_Fj), which represents the change in the total number of vesicles over time at each state. At rest, a vast majority of vesicles lay in the D state. The effect of larger numbers of molecular states on transmission was analyzed by adding states with corresponding bidirectional rate constants between the D and pP states. In the three-state model the pP state was eliminated.

The stochastic kinetic model considers that fusion requires vesicles to arrive at the P state. Since the classical kinetic differential equations do not describe correctly the collective kinetics of a small number of vesicles (∼10,000 as compared to Avogadro’s number), we used instead the master Equation 1 for the probability distribution P(x, t; x₀, t₀) (Gillespie, 1976), whose solution describes the temporal evolution of the six transition probabilities between kinetic states. The rate constants are conventional probabilities per time unit (Gillespie, 1992):

The solution of Equation 1 was simulated using the Gillespie algorithm (Gillespie, 1976), which emulates random transitions connecting different X(t) states. The fundamental equation of the Gillespie algorithm for the time evolution of the system is:

Equation 2 predicts the probability that at a state X(t) = x, the next kinetic transition R_j, will occur at the next infinitesimal time [t + τ, t + τ + dτ]. The random continuous variable τ advances the time in the simulations by the amount:

with r₁ being a random number distributed uniformly in the interval (0, 1).

The probability distribution p(j, τ) dτ mimics the solution of the stochastic kinetic Equation 1 and plays a key role in the implementation of the stochastic algorithm. Thus, the random trajectories that connect different kinetic states, X(t) = x, describe the kinetic evolution of the vesicle pool.

The algorithm for the kinetic sequence can be summarized as follows: (1) The simulation begins by setting the initial state of the system X_o at time t_o. (2) The propension functions a_j(x) and their sum a_o(x) = Σa_j(x) are calculated for each different time t. (3) The values of the discrete random variables j is chosen as the smallest integer that satisfies, ∑k=1jak(x)>r2ao with r₂ a random number distributed uniformly in the interval [0,1]. The continuous random variable τ is generated by applying Equation 3. (4) The transition to the next kinetic state x → x + v_j and the time shifts to t → t + τ are calculated. (5) A new state (x, t) is obtained, and the procedure returns to step (1).

The simulation starts with No = 10,000 vesicles accumulated in the D state. In such conditions X(t = 0) = Xo = (D(t = 0) = No, and pP(t = 0) = 0, P(t = 0) = 0, F(t = 0) = 0). As the simulation progresses, the distribution of vesicles among the different states becomes stationary in about 5 min of the simulation. After this time our measurements in the simulations are made.

The activation energies involved in the molecular transitions from docking to exocytosis lay in the same order of magnitude (Li et al., 2007; Sudhof and Rothman, 2009). Therefore, we initially considered that α₁ = α₂ = α₃ = α, and β₁ = β₂ = β. This strategy proved successful for reproducing every release mode. The α value used in cat simulations was estimated from the frequency distribution of spontaneous miniature potentials (Boyd and Martin, 1956a,b). The β and ρ values were fitted independently. Once adequate fittings were obtained, the variable values were evaluated within two logarithmic units. The model was simplified by using the coefficient λ = β/α, which permitted to evaluate the kinetic behavior in terms of the relative magnitudes of α and β. The code used in this study is available in the following repository: https://github.com/alexini-mv/kinetic-neurotransmission.

Presynaptic calcium elevations upon brief depolarization were modeled by adding a function f(t) to the forward rate constants, which acquired the form α_s = α + f(t). The kinetics of the calcium current decay in squid giant synapse experiments (Llinás et al., 1981a,b) served as the baseline. The onset of calcium transient was considered as instantaneous for the calcium channels in presynaptic neuromuscular terminals that are tightly bound to the fusion complex (Harlow et al., 2001; Nagwaney et al., 2009). Adjustments in the amplitude (in arbitrary units) and decay time (ms) of the artificial calcium elevation rendered successful results.

For our simulations it was more convenient to express the decay time τ_e of the calcium elevation instead of the decay time of the current, since according to the residual calcium hypothesis (Katz and Miledi, 1968; Kamiya and Zucker, 1994; Matveev et al., 2006), it is the residual free intracellular calcium after the impulse that promotes facilitation. The decay time of the calcium elevation was defined as:

where t_s is the stimulation time. The τ_e value was adjusted for each experimental protocol in the range of 0.05–1.5 ms. Once adjusted, the parameters of the calcium signal remained the same for each experiment. Calcium currents in certain central synapses may facilitate or depress upon subsequent stimulation (Borst and Sakmann, 1998; Cuttle et al., 1998; Forsythe et al., 1998; Inchauspe et al., 2004; Ishikawa et al., 2005; Xu and Wu, 2005; Mochida et al., 2008). However, our model rendered accurate results without any such modulation.

The simulations were made in a custom-designed code using Python 3. Calculations were carried out in a personal computer with an AMD Ryzen 5 2500U processor.

The Pearson significance was calculated using a routine of the software Mathematica: https://reference.wolfram.com/language/ref/DistributionFitTest.html.