In this paper we consider two conceptual “release-site” models (Sterratt et al., 2011) for short term synaptic depression, due to stochastic unavailability of a vesicle:

Availability Model 1: In this model it is assumed that once a specific vesicle is released, the time at which it is next available for release is a random variable that depends only on the time since it was released. This random variable is not affected by subsequent pre-synaptic spikes.

Note that these models treat a single vesicle as if it is a conserved object that switches between two states. Of course in reality the vesicle is not conserved, and a more accurate description is to say that a vesicle release site that can contain at most a single vesicle either (1) does contain a vesicle, or (2) does not contain one.

Below we show that Availability Model 1 and Availability Model 2 are mathematically equivalent, given an assumption that the random variable describing availability times is exponential. This is a standard assumption, because it provides good fits to experimental data, and therefore underpins models developed in conjunction with experimental data on short term depression (for example, Tsodyks and Markram, 1997). However, it is feasible that better fits to data might discard the exponential assumption, and in that case it would be necessary to consider how stochastic simulations need to differ for each model. As we show below for a non-exponential example (Figure 6), the two models provide significantly different outcomes.

In this paper we consider two conceptual models for the stochastic release of a single vesicle upon arrival of a pre-synaptic spike:

Release Model 1: In this model it is assumed that if the single vesicle is available, then it is released with a constant probability upon arrival of a pre-synaptic spike, and this probability does not change over time.

Release Model 1 is a classical model of probabilistic release (Vere-Jones, 1966). Release Model 2 is appropriate when a synapse is known to exhibit facilitation. Usually, based on experimental evidence (see, for example, Markram et al., 1998), the change in release probability (given availability) over time is modeled as increasing by a percentage of the current probability of non-release, whenever a pre-synaptic spike arrives (usually independently of whether a vesicle is available or released) and then decaying exponentially over time to a constant rest probability, for as long as no more pre-synaptic APs arrive (Tsodyks et al., 1998).

Release Model 2 is also appropriate when a synapse is known to exhibit a different form of short term depression to that modeled by the lack of vesicle availability. In this type of depression, known as “release-independent depression,” the probability of vesicle release (given its availability) is reduced by arriving pre-synaptic spikes independently of whether the vesicle is released, due to different mechanisms from those that cause facilitation (Fuhrmann et al., 2004). In some models, facilitation and release-independent depression are assumed to be present simultaneously (Graham and Stricker, 2008; Scott et al., 2012).

A single vesicle obviously cannot be released if it is not available, but it is assumed that an available vesicle remains available until released. This is the key feature of the conceptual models we study where both availability and release are modeled as stochastic.

In this paper, when we consider a conceptual model where there is a pool consisting of at most N vesicle release sites, each containing at most a single vesicle, we use the typical assumption that the release and availability of each single vesicle occurs independently of that in the other vesicle release sites. Note that although this model is typical, it may not always be accurate (Quastel, 1997).

In this paper, when we consider a conceptual model where there are N repeated trials for the same sequence of pre-synaptic APs, and a single vesicle in a single release site, we assume that the availability or release of the vesicle is independent for each trial.

Note that the outcome for a model where there are N such repeated trials is equivalent mathematically to a conceptual model where there is a pool of N vesicle release sites with at most a single vesicle available, for a single trial of the sequence of APs.

Since an experimental protocol is more amenable to studying repeated trials for a single release site and the same sequence of APs, we will refer to the case of N trials rather than N release sites.

The content of this paper regarding stochastic simulations can be extended to a scenario where multiple vesicles are available in a release site, and also where multiple such sites are available, potentially each with different numbers of vesicles. However, we do not discuss this further, as the most important observations are relevant to sites containing single vesicles. Further discussion of evidence for multiple release sites can be found in Loebel et al. (2009).

There is a specific random variable that exists in both conceptual availability models: the time taken for vesicle to become available following a successful release at time t = t_s. We label this random variable as T_a1 for Availability Model 1, and as T_a2 Availability Model 2.

There is a specific random variable that exists in both conceptual release models: the event that a vesicle is released, or not released, upon arrival of the i–th pre-synaptic AP at time t = t_{AP, i}. We label this random variable as R(t_{AP, i)}. This random variable is binary, it exists only at each AP time, it depends on the last time at which a vesicle was released, t_s, and we denote its outcomes as α if a vesicle is released and as β if it is not. We denote the probability that the event α occurs at time t, given the vesicle is available, as P_r|a(t). The random variable has a probability mass function, and this is given for Availability Model 1 by (6)Prob(R(tAP,i)=α|ts)=Pr|a(tAP,i)FTa(tAP,i−ts);Prob(R(tAP,i)=β|ts)=1−Pr|a(tAP,i)FTa(tAP,i−ts), where t_{AP, i} > t_s, and for Availability Model 2 by (7)Prob(R(tAP,i)=α|ts)=Pr|a(tAP,i)Pa,2(tAP,i|ts);Prob(R(tAP,i)=β|ts)=1−Pr|a(tAP,i)Pa,2(tAP,i|ts). Note that in Release Model 1, P_r|a has no time dependence [i.e., P_r|a(t_{AP, i}) = P_r|a], but this is the only difference in comparison with Release Model 2 (see Equation 1). Consequently, provided the release probability has been calculated correctly at each point in time during a simulation, there are no other differences in a stochastic simulation implementation in comparison with Release Model 2.

A relevant binary-valued stochastic process can be stated based on the random variables described above, namely, the process describing whether a vesicle is available at any point in time. A succinct description of this process is given in Loebel et al., (2009, Equation 4), where it is expressed in terms of a differential equation. Following the notation used in that description, we label the stochastic process as σ(t), and let σ(t) = 1 when the vesicle is available and let σ(t) = 0 otherwise. The process is fully described by the following equation: (8)dσ(t)dt=−σ(t−)R(t)δ(t−tAP,i)             +(1−σ(t−))δ(t−tAP,i−Ta), where the notation ⁻ in t⁻ is used as shorthand to represent t⁻ = t − ∈, where ∈ is a very short time period; thus when t = t_{AP, i} then t⁻ is the time instant immediately prior to AP i arriving. We have assigned α = 1 and β = 0 as the possible values of the random variable R(t). During intervals of time for which σ(t) = 0, the right hand side of Equation (8) is just δ(t − t_{AP, i} − T_a), and mathematically, the remaining terms in Equation (8) describe the fact that σ(t) can jump from 0 to 1 only at the time t = t_{AP, i} + T_a. Similarly, Equation (8) is such that σ(t) can jump from 1 to 0 only when both R(t) = α = 1 and t = t_{AP, i}, or equivalently, R(t_{AP, i}) = 1, which means a vesicle is released when AP i arrives.

Note that in Loebel et al. (2009), the event where σ(t) jumps from 0 to 1 is stated to be modeled as a Poisson process. A Poisson process has exponentially distributed times between events, and therefore the conceptual model in Loebel et al. (2009) is in this sense the same as our Availability Model 1 with exponentially distributed T_a1, with mean τ_a. However, for an actual Poisson process, events will continue to occur for all time, not just when vesicles are currently unavailable, which is at odds with our stated conceptual model. Nevertheless, it can be inferred that in Loebel et al. (2009) that Poisson events are ignored when σ(t) = 1.

Does this mean that the distribution of times until a vesicle becomes available is different in each conceptual model, since in the Poisson process, the exponential time to arrival begins at the time of the previous Poisson event, whereas in Availability Model 1 begins at the most recent release time? The answer is no, due to the independence of events in Poisson processes (the same reason that Availability Models 1 and 2 are equivalent for exponentially distributed availability times). Therefore, there will be no difference when a stochastic simulation implementation of the Loebel et al. (2009) conceptual model is carried out, compared with an implementation of our Availability Model 1 with exponentially distributed arrival times. However, if the arrival times are not exponential, and the corresponding non-Poisson process replaces the Poisson process in the Loebel et al. (2009) conceptual model, the results will not be the same.

Differential equation notation is often used to express how the mean fraction of available vesicles, N_a(t), changes over time in two ways: either upon a spike arrival, or between spike arrivals (Tsodyks and Markram, 1997; Fuhrmann et al., 2002; Scott et al., 2012). The typical form of such expressions is dNa(t)dt=1−Na(t)τa−Nr|aNa(t−)∑i=1Kδ(t−tAP,i), where t_{AP, i} is the arrival time of the i–th AP, out of a total of K, and N_r|a is the mean fraction of available vesicles released by the i–th AP. This differential equation can be easily solved in closed form (e.g., Tsodyks and Markram, 1997) to get (9)     Na(t)=1,                             t<tAP,1,Na(tAP,i)=(1−Nr|a)Na(tAP,i−),   t=tAP,i,     Na(t)=1−(1−Na(tAP,i))     t∈[tAP,i,tAP,i+1),                 exp(−(t−tAP,i)/τa) where i = 1, …, K.

Note that the change over time in the fraction of trials in which the vesicle is available clearly has a dependence on both (1) the time since the most recent pre-synaptic AP and (2) on the fraction of vesicles available at the time of the most-recent pre-synaptic AP. Consequently, the deterministic mean model should be interpreted as explicitly solving for the conditional mean number of vesicles released at each AP arrival, given the number that are available for release.

Remark 1: It is clear that the mean model accurately reflects Availability Model 2 generally, and in the specific case stated above, assumes exponential availability times following each AP arrival. Moreover, we have discussed that Availability Models 1 and 2 are equivalent for exponentially distributed availability times, and hence the stated mean model also accurately reflects Availability Model 1 for this specific case.

Remark 2: The deterministic mean model does not, however, accurately reflect Availability Model 1 for non-exponentially distributed availability times, since under Availability Model 1, the fraction of trials in which a vesicle should be released, given that it is available, should be based on the trial-dependent time since a vesicle was released, not solely on the time since the most recent AP. Therefore, the right hand side of a differential equation describing the mean number of trials in which a vesicle is available should have an additional term for each AP that occurs prior to the current AP. Moreover, if each additional term describes the mean number of trials in which vesicles have not become available since the i–th AP, the results will potentially become increasingly inaccurate with the time elapsed since the i–th AP.

One possibly useful element in any extension of mathematical analysis to this case of multiple trials might be an iterative expression articulated in a different context by McDonnell et al. (2002, 2008), that can be adapted to describe the conditional probability that u vesicles are available across Z trials, even if the time they were released differs. This approach does not suggest a straightforward method for implementing a stochastic simulation, but as described by McDonnell et al. (2002, 2008), there are simple expressions for the conditional mean and variance, and these could potentially be used within a deterministic equation that describes how the mean number of trials in which a vesicle is available changes with time.

In section 5, we compare the results of stochastic simulations with results for the mean obtained from Equation (9). We also use a result for a scenario where pre-synaptic APs arrive at the synapse periodically with frequency f Hz so that the AP times are t_{AP, i} = i/f, i = 1, 2, …. In this case, it is well known that the mean fraction of vesicles available quickly decays to a constant steady state value, N⁻_ss := N_a(t⁻_{AP, i + 1}) = N_a(t⁻_{AP, i}). As shown in (Abbott et al., 1997; Matveev and Wang, 2000b), this can be obtained from Equation (9) (which hold for Availability Model 2 generally, and for Availability Model 1 with exponentially distributed availability times) to get the mean fraction of vesicles available for release just prior to a pre-synaptic spike as (10)Nss−=(1−exp(−1/(fτa)))1−exp(−1/(fτa))(1−Nr|a).

The following pseudo-code illustrates how simulations of the random variables described above can be implemented in stochastic simulations.

Correct Implementation 1, for AM1 Set: NextAvailabilityTime = 0 For each pre-synaptic spike, i=1:K, occurring at time t_i if t_i >= NextAvailabilityTime //vesicle is available if Pr_given_a(t_i) > unifrand() //Release the vesicle Set: LastReleaseTime = t_i //reset the next availability time, for //exponentially distributed availability times NextAvailabilityTime = LastReleaseTime +exprand(tau) end end end //unifrand() generates a uniformly distributed random number between 0 and 1 //exprand() generates an exponentially distributed random number with mean tau

The pseudo-code variable LastReleaseTime represents our mathematical variable, t_s. A direct translation of this pseudo-code into the probability that the vesicle will be released upon the arrival of AP i, given t_s, obtains P_r|a(t_{AP, i})Prob(t_{AP, i} > t_s + T_a). This can be expressed as P_r|a(t_{AP, i})Prob(T_a < t_{AP, i} − t_s) = P_r|a(t_{AP, i})F_Ta(t_{AP, i} − t_s), and thus exactly matches Equation (6), as required.

There are also several ways in which the conceptual model has been, or could be, erroneously translated into a stochastic simulation, and these are described in the following subsections.

The following pseudo-code illustrates how simulations of Availability Model 2 can be implemented in stochastic simulations. Note that the only difference in comparison with the pseudo-code for Availability Model 1 is that the time of next availability (for unavailable vesicles only) is dependent only on the last spike arrival time, not the last release time, in order to match the conceptual model.

Correct Implementation 1 for AM2 Set: NextAvailabilityTime = 0 For each pre-synaptic spike, i=1:K, occurring at time t_i if t_i >= NextAvailabilityTime //vesicle is available if Pr_given_a(t_i) > unifrand() //Release the vesicle and reset the next availability time, for // exponentially distributed availability times NextAvailabilityTime = t_i +exprand(tau) end else //vesicle is unavailable; reset the next availability time, for //exponentially distributed availability times NextAvailabilityTime = t_i +exprand(tau) end end

A direct translation of this pseudo-code into the probability that the vesicle will be released upon the arrival of AP i, given that it was not released by the time of AP i − 1, obtains P_r|a(t_{AP, i})Prob(t_{AP, i} > t_{AP, i − 1} + T_a). This can be expressed as P_r|a(t_{AP, i})Prob(T_a < t_{AP, i} − t_{AP, i − 1}), and thus exactly matches Equation (7), as required, upon substitution of Equation (5).

It is possible to incorrectly implement Availability Model 2 in a manner directly analogous to that in the first incorrect implementation of Availability Model 1. However, an implementation analogous to the second incorrect implementation of Availability Model 1, will actually be correct for Availability Model 2, since now the probability of availability is dependent only on the last AP arrival time.

We have aimed in the pseudo-code implementations above to describe computationally efficient algorithms that require as few random numbers to be generated as possible.

We note that the implementation suggested, for example, in Loebel et al. (2009) [see also Sterratt et al., (2011, p. 188)] involves an accurate approximation of a true Poisson point process, and this approximation is particularly relevant to any simulation in which time is discretised into uniform intervals of Δt, such as in most simulations that involve numerical solution of differential equations. The well-known approximation states that provided that Δt « τ_a, a Poisson point process event occurs within any given time interval of duration Δt with probability Δtτa.

A stochastic simulation based on this approximation requires comparison of Δtτa with a uniform random number at every time step of the simulation between times 0 and t_K. It is possible to alter the implementation so that the Poisson events are only calculated during the simulation, rather than prior, where a comparison of a uniform random number with Δtτa is carried out for every time step following vesicle release, until a random number is generated that is larger than Δtτa.

However, such implementations are potentially very inefficient, because many random numbers must usually be generated for every unavailable vesicle, whereas only one random number need be generated in, for example, Correct Implementation 1 for AM1.

We consider a scenario where pre-synaptic APs arrive at a synapse periodically with frequency f Hz. We consider Z repeated trials following an initial condition where a vesicle is assumed to have just been released, in all trials, at the start of our simulations. We calculate the number of trials in which the next vesicle release occurs after the first spike, the second spike and so forth. We obtain results for f between 5 and 150 Hz, for Z = 100,000, τ_a = 0.5 s and K = 50 pre-synaptic spikes (as a maximum; the simulation stops when a vesicle is first released). Thus, the AP times are t_{AP, i} = i/f, i = 1, 2, …,50.

We estimated the probability that the vesicle was next released after i = 1, 2, … 20 APs following vesicle release at time t = 0, by evaluating the fraction of trials in which the vesicle was first released after the i–th spike. We then calculated the absolute value of the difference in the estimated probability for several correct and incorrect implementations, and also the relative error, relative to the correct version.

In order to clarify the significance of the values we obtained for absolute and relative error, we also considered a simulation where H = 100 spikes per trial periodically arrive with frequency f, and for each implementation counted the total number of vesicles released as a function of f. We then compared the mean number released after H spikes, calculated from Z = 10,000 repeats of each implementation, as well as the maximum and minimum numbers released.

Finally, in order to show that our simulations and mathematical analysis is correct for Availability Model 1 for both periodic and non-periodic AP arrivals, we compare correct and incorrect implementations for each case with results predicted by the Equations (9) and (10).

In order to demonstrate that a non-exponential availability model provides different outcomes for Availability Models 1 and 2, we consider the case of Rayleigh distributed availability times, with mean τ_a = 0.5 s.

Figure 6 shows the fraction of 1000 trials in which vesicles are released, for both availability models, in response to a sequence of 20 periodically arriving APs, with frequency 10 Hz, and to a sequence of 50 APs arriving at times corresponding to a Poisson point process, with mean frequency 10 Hz.

The results shown in Figure 6 were obtained both by a direct adaptation of the correct stated pseudo-code above to Rayleigh distributed availability times, and also by direct adaption of the binomial approaches described, for Availability Model 2. The binomial results for Availability Model 1 required a more complex algorithm, where the number of vesicles unavailable due to release from all previous spikes needed to be tracked, and consequently many more binomial random numbers generated than for Availability Model 1. The data can be seen to match in either implementation, but to be quite different for each Availability Model.

We have shown that various correct implementations of a stochastic simulation of either Availability Model 1 or 2 are possible. However, it is also possible to incorrectly implement either model. For Availability Model 1, two kinds of incorrect implementation result in more vesicle releases than should be the case. For Availability Model 2, an incorrect implementation results in less vesicle releases than should be the case.

We have also shown that some correct implementations are more efficient than others. In particular, we first stated an implementation that requires only a single random number to be generated each time a vesicle is released. This is more efficient than an implementation based on generation of a Poisson process that determines availability times, and much more efficient than generating a random number for every time step in a simulation.

We also have shown that when multiple independent vesicle releases are considered, the most efficient stochastic simulation implementation can be achieved by generating binomial random numbers.

We have discussed that the two availability conceptual models are equivalent when the availability times are exponentially distributed, and this can be derived as a consequence of a well known property of a homogeneous Poisson point process. When the availability times are non-exponential, the two availability conceptual models we consider are generally non-equivalent.

As we have shown, these points are important in terms of their consequences for the implementations that can be used to correctly simulate Availability Model 1. Another consequence is that the popular differential equation approach to describing the mean number of available vesicles could have analogous correct simple forms for Availability Model 2, but not for Availability Model 1.

When T_a is not exponentially distributed, no simple adaptation of the binomial approach will work for Availability Model 1, because there is no simple procedure for calculating random numbers corresponding to the random variable W. However, it is possible to keep track of how many vesicles were made unavailable by each AP, and how many are restored by the time of each subsequent AP, and perform a stochastic simulation that calculates as many binomial random numbers as there have been prior APs for which unavailable vesicles still exist. This procedure must make use of Equation (2) to calculate the probability of availability by the time of the next spike, using the time of the previous spike.

Such an algorithm may be more efficient than independently simulating each trial, provided that APs arrive relatively slowly compared with the mean time for a released vesicle to become available. Figure 6 shows that an implementation of this approach for Availability Model 1 and Rayleigh distributed available times agrees with data from an approach that individually simulates each trial.

For Availability Model 2, the binomial approach will work for arbitrary F_Ta, since the probability of an unavailable vesicle becoming available will be the same for all trials, as was the case for the data in Figure 6. Moreover, only the cumulative distribution function of the availability times need be known to carry this out, whereas Equation (2) needs to be computed for Availability Model 2. Such an algorithm is a straightforward extension of an algorithm described by Quastel (1997) for constant probabilities for any released vesicle to be available by the next spike.

We have considered only the simplest models in this paper. Other more complex models of availability and release have been proposed. For example, the mean times to availability for a vesicle may change over time (Wong et al., 2003), vesicles may also be released spontaneously in the absence of pre-synaptic APs (Sterratt et al., 2011), and multiple vesicles may be readily available for release at any site (although in this model only one can be released per pre-synaptic spike) (de la Rocha and Parga, 2005). Stochastic simulations that faithfully reflect these models can be readily devised by extension of the algorithms presented in this paper.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.