Bayesian Inference of Synaptic Quantal Parameters from Correlated Vesicle Release

1. Introduction

The statistics and dynamics of stochastic synaptic filtering determine how information is communicated between neurons. Synapses act as activity-dependent filters on the transfer of neuronal signals, suppressing or amplifying trains of inputs to the postsynaptic cell relative to isolated stimuli, in a phenomenon known as short-term plasticity or synaptic dynamics (Zucker and Regehr, 2002; Abbott and Regehr, 2004; Mongillo et al., 2008). An action potential in the presynaptic cell triggers an influx of Ca²⁺ into synaptic terminals, causing a probabilistic all-or-none release of neurotransmitter at each active vesicle docking site on the presynaptic membrane. The neurotransmitter binds to channels on the postsynaptic cell resulting in, for example, an excitatory post-synaptic potential (EPSP) “built up statistically of the all-or-none events that are similar in size and distribution to spontaneous miniature” postsynaptic potentials (del Castillo and Katz, 1954). Depletion of vesicles available at active sites can cause an activity-dependent reduction in synaptic efficacy (Eccles et al., 1941) whereas a build-up of Ca²⁺ in the presynaptic terminal can increase the probability of neurotransmitter release (Dudel and Kuffler, 1961). Synaptic transmission is thus both fundamentally stochastic (del Castillo and Katz, 1954; Fatt and Katz, 1952; Stein, 1965) and history dependent (Furukawa et al., 1982; Abbott, 1997; Tsodyks and Markram, 1997).

Initial analyses of paired-cell data used the amplitude distribution of isolated EPSPs to identify quantal peaks corresponding to sums of similar mini amplitudes (Boyd and Martin, 1956; Liley, 1956; Kuno, 1964; Kuno and Weakly, 1972; Bennett and Florin, 1974; Bekkers, 1994); for a review see Bennett and Kearns (2000). While this was an effective approach for extracting the properties of neuromuscular synapses (del Castillo and Katz, 1954) the greater variation in mini amplitudes at central synapses (Hanse and Gustafsson, 2001; Franks et al., 2003; Hardingham et al., 2010) necessitated different techniques to recover robust results in the central nervous system. Mean-variance analysis was developed to obtain estimates of the maximum number of vesicles that can be released by a single stimulus (Silver et al., 1998; Clements, 2003; Silver, 2003). Initial applications relied on conducting experiments under a variety of conditions, in particular varying the extracellular Ca²⁺ concentration to alter the vesicle release probability (Foster and Regehr, 2004; Birò et al., 2005). Brémaud et al. (2007) and Loebel et al. (2009) increased the practicality of the method by using short-term vesicle depletion to vary the effective release probability under a single experimental condition. Their analyses showed that multiquantal release underlies the wide range of EPSP amplitudes observed (Song et al., 2005; Lefort et al., 2009) and that, in general, it is not the case that the number of distinct anatomical contacts equals the maximum number of readily-releasable vesicles as was put forward by the single-vesicle hypothesis (Kuno, 1971; Korn et al., 1981).

More recent approaches have introduced a principled Bayesian approach to infer synaptic parameters. Bayesian inference determines the extent to which experimental evidence supports a given set of model parameters. This relies on the fact that the probability of a certain model being correct given observed data is proportional to the probability of observing that data given that the model is correct. As such it makes maximal use of data, including every observation rather than extracting moments as in previous approaches. This framework was first applied to neurophysiological synaptic data by Turner and West (1993) to extract the number of components in a unitary EPSP. More recently, McGuinness et al. (2010) used Bayesian analysis to measure presynaptic Ca²⁺ concentrations and Bhumbra and Beato (2013) used an exact Bayesian approach to extract quantal parameters from measurements of isolated EPSPs, demonstrating that accurate parameter estimates could be obtained from less data than with existing mean-variance methods.

Inference on isolated EPSPs, however, does not allow recovery of synaptic parameters associated with short-term plasticity. Costa et al. (2013) addressed this issue in a Bayesian framework using the Tsodyks-Markram model of short-term plasticity (Tsodyks et al., 1998) with a likelihood that approximated synaptic amplitude distributions during patterned input as uncorrelated Gaussians around the mean amplitudes. Though this approach does not account for correlations between closely-timed synaptic events, the method nevertheless allowed for accurate inference of a number of synaptic parameters. However, correlations between successive PSPs, which can be significant even at stimulation rates below 10 Hz, (del Castillo and Katz, 1954; Thomson et al., 1993; Fuhrmann et al., 2002) can provide a useful source of additional information for inferring model parameters. This is particularly the case for data sets that feature only a few repreated stimulations or only one series of patterned PSPs such as would be the case for spontaneous in-vivo recordings.

The main barrier to extending the Bayesian approach to a model that allows simultaneous recovery of both quantal and dynamic properties is the calculation of the likelihood of seeing a particular train of amplitudes in response to a certain pattern of presynaptic stimuli. This probability is dependent on the correlated vesicle releases during previous events and the number of possibilities therefore grows exponentially with the number of PSPs. Naively, this would appear to make the problem intractable. However, two independent studies (Barri et al., 2016; Bird, 2016) recently provided a solution to this problem by exploiting the underlying Markovian nature of the problem thereby allowing for the computation of the exact probability of a given set of observed amplitudes with a complexity that grows only linearly with PSP number. Here we develop the method, originally presented in Bird (2016), to show how the likelihood may be written in a compact mathematical form as a matrix product. This allows for efficient calculation of the posterior distribution from which, for example, the covariance of the inferred parameters can be analyzed. Our complete Bayesian method may be thought of as combining the method for inferring quantal parameters for isolated PSPs developed by Bhumbra and Beato (2013) with the method for inferring mean synaptic dynamics (without including correlations) developed by Costa et al. (2013). As well as describing the mathematical solution we additionally provide the software code to perform Bayesian inference for a variety of models of synaptic dynamics as part of this publication.

2. Methods

In this section we define the general class of synaptic models our inference procedure applies to before specifying a commonly used depression-facilitation model of neurotransmitter release that will be used for illustrative purposes. The coupling of the presynaptic model to the post-synaptic voltage response is then defined.

2.1. The Class of Synaptic Dynamics Models

The method presented here is applicable to a broad class of synaptic models. The synapses this method can be applied to are assumed to have a number n of vesicle release sites to which neurotransmitter vesicles can dock. On arrival of the mth presynaptic spike at time t_m neurotransmitter is released independently from each docked vesicle with probability u_m. The binary occupancy variable x(t) for single release site obeys

$\begin{array}{ll}\frac{dx}{dt}=(1-x){\displaystyle \sum _{\left\{t_r\right\}}}\delta (t-t_r)-x{\displaystyle \sum _{\left\{t_m\right\}}}{\delta }_m\delta (t-t_m)& (1)\end{array}$

where t_r are restock events (which occur at a rate that may be dependent on the presynaptic action potential times) and δ_m is a binary random variable signifying release of neurotransmitter that is equal to 1 with probability u_m and 0 otherwise. The stochasticity in t_r and δ_m is considered to be statistically independent across the n vesicle release sites. Note also that in this formulation any dynamic quantity (such as x(t)) multiplying a Dirac-delta function is evaluated just before the arrival of the impulse. The expected change in occupancy before and after a presynaptic action potential can be straightforwardly derived to give

$\begin{array}{ll}{\langle x\rangle }_m^{\oplus }={\langle x\rangle }_m^{\ominus }-u_m{\langle x\rangle }_m^{\ominus }& (2)\end{array}$

where ${\langle x\rangle }_m^{\ominus }$ is the probability that a release site is occupied just before and ${\langle x\rangle }_m^{\oplus }$ just after the mth spike. Similarly, the probability of occupancy just before (m + 1)th AP can be related to the occupancy just after the mth AP as

$\begin{array}{ll}{\langle x\rangle }_{m+1}^{\ominus }=1-(1-{\langle x\rangle }_m^{\oplus })(1-g_m)& (3)\end{array}$

where g_m is the restock probability. For certain models g_m can depend on the history of the presynaptic APs. Together the recursion relations (3) and (5) give the occupancy probability for an arbitrary train of presynaptic action potentials. The initial condition is typically taken as being ${\langle x\rangle }_1^{\ominus }=1$, where all release sites are stocked. These dynamics cover a range of models such as vesicle depression (Tsodyks and Markram, 1997), depression with facilitation (Varela et al., 1997; Tsodyks et al., 1998; Fuhrmann et al., 2002), frequency-dependent recovery (Fuhrmann et al., 2004) and augmented recovery (Wang and Kaczmarek, 1998; Hosoi et al., 2007). For an in-depth discussion, see Appendix A.

2.2. Illustrative Synaptic Model with Depression and Facilitation

To provide an example of the method we use a commonly used model that combines a depression mechanism caused by vesicle release and a constant restock rate with a facilitation mechanism that models the effect of increased release probability due to transient increases in calcium concentrations in the presynaptic terminal (Varela et al., 1997; Tsodyks et al., 1998; Fuhrmann et al., 2002). The restock process is Poissonian and has constant rate 1/τ_D, where τ_D is commonly referred to as the depression time constant; therefore the restock probability required for Equation (3) is simply

$\begin{array}{ll}{g}_m=1-e^{-T_m/{\tau }_D}& (4)\end{array}$

where T_m = t_m+1 − t_m is the time between the mth and (m + 1)th APs. Let p₀ be the baseline value of the probability of release, and p₁ be the facilitated release probability immediately after an isolated spike. Let u(t) be the time-dependent release probability. In the absence of stimulus, u(t) decays back to p₀ with timescale τ_F. The dynamics of u(t) therefore obeys

$\begin{array}{ll}\frac{du}{dt}=\frac{1}{{\tau }_F}(p_0-u)+(1-u)\left(\frac{p_1-p_0}{1-p_0}\right){\displaystyle \sum _{t_m}}\delta (t-t_m)& (5)\end{array}$

where the (1 − u) prefactor of the Delta functions prevents the probability going above unity. In this setup u = p₀ if the previous spike was a long time ago, then on the arrival of a spike it jumps to u = p₁. Because it is a facilitation model we have p₀ < p₁ < 1. Note that this formulation of parameters allows the facilitated release probability p₁ to be fixed independently of the initial release probability p₀ and maps directly to the original quantal facilitating and depressing synaptic model of Fuhrmann et al. (2002) with p₀ = U_SE and p₁ = U_SE + (1 − U_SE)U₁ using that paper's notation. The values of u(t) just after the mth and before the (m + 1)th action potentials ($u_m^{\oplus }$ and u_m+1 respectively) are defined by the following recursion relations

$\begin{array}{l}{u}_m^{\oplus }=u_m+(1-u_m)\left(\frac{p_1-p_0}{1-p_0}\right)\mathrm{\text{ and}}\\ u_{m+1}=p_0+(u_m^{\oplus }-p_0)e^{-T_m/{\tau }_F}& (6)\end{array}$

where the initial conditions are that u₁ = p₀ and $u_1^{\oplus }=p_1$. This gives the release probability before each presynaptic spike required for Equation (2). The dynamics of the restock probability g are unaffected and are given by Equation (4). A special case of this model that has one less free parameter is when the release probability doubles after an isolated spike and so p₁ = 2p₀ (Tsodyks et al., 1998).

2.3. EPSP Amplitude Distribution

The post-synaptic amplitude statistics for single vesicle release of neurotransmitter is modeled by a gamma distribution with mean μ_a and standard deviation σ_a (μ_a is assumed to be greater than σ_a to ensure that a zero amplitude is not the most likely). This is preferred over a normal distribution on empirical grounds and ensures that amplitudes are always positive (Robinson, 1976; Hanse and Gustafsson, 2001; Bhumbra and Beato, 2013). However, it is reasonable to assume that background noise is normal with standard deviation σ_b and is independent of EPSP amplitude. Note that this choice of amplitude generation is identical to that described for isolated EPSPs in (Bhumbra and Beato, 2013). With this choice, if k vesicles release neurotransmitter from among the n possible release sites, the observed EPSP amplitude A is written A = ψ + ϕ where ψ is the release-dependent component and ϕ the independent Gaussian noise. Because ψ is the sum of k individual quantal amplitudes, each of which are independently identically gamma distributed, its distribution is also gamma-distributed with

$\begin{array}{llll}\text{P}\left[\psi \right]=\frac{{\lambda }^{\beta }}{\Gamma (\beta )}{e}^{-\lambda \psi }{\psi }^{\beta -1}& \mathrm{\text{where}}& \beta =k\frac{{\mu }_a^2}{{\sigma }_a^2}\text{ and }\lambda =\frac{{\mu }_a}{{\sigma }_a^2}.& (7)\end{array}$

The distribution for the measured EPSP amplitude A, given k release events, is therefore a convolution between the gamma and normally distributed components of the noise

$\begin{array}{ll}\text{P}\left[A|k\right]=\frac{{\lambda }^{\beta }}{\Gamma (\beta )}\frac{1}{{(2\pi {\sigma }_b^2)}^{\frac{1}{2}}}{\displaystyle {\int }_0^{\infty }}dy{e}^{-\lambda y}{y}^{\beta -1}{e}^{-\frac{{(A-y)}^2}{2{\sigma }_b^2}}.& (8)\end{array}$

An approach for numerically calculating this integral efficiently is provided in Appendix B.

2.4. Computational Methods and Code

An exhaustive grid-based derivation of the likelihood function for the depression-only model (see Appendix) is just within the capabilities of easily accessible computers at the time of writing. However, for more involved models with a greater number of parameters this becomes impracticable and a Markov Chain Monte Carlo (MCMC) approach was used instead. Here priors are taken to be flat (uninformative) for all parameters for illustrative reasons: more informative priors can be included as required. For the MCMC implementation, parameter space is discretised into a grid and the sampler is initialized at a random point consistent with any restrictions on the model parameters. Moves are proposed to each adjacent grid point with equal probability and accepted or rejected based on the log-likelihood ratio of the current and proposed points. Convergence of the sampler was examined by comparing the distributions resulting from chains initiated in different locations. It is straightforward to extend this transparent implementation in our code to include more sophisticated methods such as slice sampling. We provide MATLAB and JULIA code for the Bayesian inference of synaptic parameters from measurements of synaptic amplitudes using the Metropolis-Hastings sampling method (Metropolis et al., 1953; Hastings, 1970) described above as part of the Supplementary Material. The code covers the major synaptic dynamics models including: depression only, depression-facilitation, release-independent depression and frequency-dependent recovery. The models are described in the Appendix.

2.5. Synthetic and Experimental Data

To test the model we used both artificial and experimental data sets. Synthetic data with known parameters was generated from the synaptic-dynamics models and consisted of a series of stimulation times and stochastically determined EPSP amplitudes. For experimental data sets the data analyzed consisted of EPSP amplitudes combined with their arrival times. The data, comprising paired whole-cell patch-clamp recordings of layer-5 pyramidal neurons, was taken from a previous study (Kerr et al., 2013). Here data obtained in control conditions and in the presence of 100μM bath-applied adenosine was used. Presynaptic cells were stimulated with square-pulse currents of 5ms duration and magnitude sufficient to reliably induce a single action potential without causing bursting. Stimulation consisted of 10 spikes at 20 − 50 Hz with 10 s between traces ensuring sufficient time for full recovery and statistical independence for the next sweep. For each presentation of the same presynaptic stimulus the amplitudes of the individual EPSPs were extracted from the postsynaptic voltage trace using the voltage deconvolution method (Richardson and Silberberg, 2008) providing a vector of 10 EPSP amplitudes.

3. Results

In this section we first summarize the broad class of synaptic models our methodology applies to. We then describe the nature of the computational problem involved in calculating exact correlated likelihoods. We go on to show how the probability of observing a set of numbered release events for a chain of presynaptic action potentials can be calculated using a Markovian property. By coupling this result to the miniature PSP distribution, the full likelihood for an observed PSP amplitude train is then derived in the form of a matrix multiplication. Finally, we demonstrate the method on both synthetic and experimental data, recovering the shift in synaptic dynamics caused by the neuromodulator adenosine.

3.1. Synaptic Models

We consider synaptic models that are quantal, stochastic and exhibit short-term plasticity. The synaptic-dynamics models feature n sites where a vesicle can be present for release. If a vesicle is present just before the mth pulse then it is released with probability u_m. Between the mth and (m + 1)th pulses an empty vesicle site can be restocked with probability g_m. Both release (given a presynaptic AP) and restock events are independent between release sites. The probabilities themselves are deterministic in that they depend on the model parameters only and can be calculated in advance if the times of the action-potentials t_m are known. This formulation encompasses a very broad range of models of short-term plasticity.

When a vesicle is released, the size of the mini PSP it produces in the postsynaptic cell is modeled by a gamma-distributed random variable (see Methods). The mini PSPs induced by different vesicles are assumed to be independently identically distributed. The mean quantal amplitude is μ_a and the standard deviation is σ_a. In addition there is a normally-distributed background noise with standard deviation σ_b that is uncorrelated with EPSP amplitude.

For illustrative purposes we focus on a model of synaptic dynamics that features depression and facilitation (Tsodyks and Markram, 1997; Fuhrmann et al., 2002), though other models for which computer code is also provided are described in Appendix A. Activity reduces synaptic efficacy through vesicle depletion; however, the build-up of Ca²⁺ in a presynaptic terminal means that the probability of release given that a vesicle is present u is increased by presynaptic activity. Thus, the response to sustained activity can involve larger individual PSPs than the response to isolated spikes. Here, the model has a probability p₀ of release for an isolated pulse; immediately after an isolated presynaptic action potential the release probability increases to p₁. The release probability u returns to its initial value p₀ with a timescale τ_F. Empty release sites are restocked on a timescale of τ_D. The model is fully defined in Methods and its parameters are summarized in Table 1.

TABLE 1

Table 1. Table of inferred parameters (top) and dynamic variables (bottom) used in the model of synaptic dynamics.

3.2. The Nature of the Computational Problem

We now discuss the aim of Bayesian inference and the difficulties correlations cause in calculating the necessary quantities. We consider that the data is in the form of a set of presynaptic action-potential times t₁, t₂, ⋯ , t_M and post-synaptic amplitudes A₁, A₂ ⋯  A_M. The aim of the inference procedure is to calculate the probability densities of the parameters of the model θ = {n, p₀, p₁, τ_D⋯ } given the observed presynaptic action potential times {t₁, ⋯ , t_M} and postsynaptic amplitudes {A₁, ⋯A_M}. Bayesian inference utilizes the fact that the probability of a particular set of parameters being true, given some observed data, is proportional to the probability of observing that data given that those parameters are correct:

$\begin{array}{ll}\text{P}(\theta |A_M,A_{M-1},\cdots ,A_1)\propto {L}(A_M,A_{M-1},\cdots ,A_1|\theta ).& (9)\end{array}$

The term ${L}$ on the right-hand side is referred to as the likelihood function. A-priori calculating the likelihood appears computationally infeasible as naively it might be expected to grow exponentially with the number of observed amplitudes M. For example, consider a case with n possible release sites and a pair (M = 2 of presynaptic spikes. Then the likelihood ${L}$ is given by

$\begin{array}{ll}{L}(A_2,A_1|\theta )={\displaystyle \sum _{k_2=0}^n}{\displaystyle \sum _{k_1=0}^n}\text{P}\left[A_2|k_2\right]\text{P}\left[A_1|k_1\right]\text{P}\left[k_2,k_1\right]& (10)\end{array}$

where k_m is the number of vesicles released by the mth spike. Because of the nested sums there are (n + 1)² additive terms in this expansion, and more generally the number of terms in the expansion grows exponentially with the number of presynaptic action potentials ~ (n + 1)^M. Written in this form it is clear that the problem becomes computationally prohibitive for long trains of presynaptic spikes and this is what makes calculation of the likelihood difficult for the complete model. The complexity arises from the quantal part of the likelihood P[k₂, k₁]; the individual amplitudes A_m are dependent only on the number of vesicles k_m released by each action potential.

Note that if correlations are ignored and the approximation P(k₂, k₁) ≃ P(k₂)P(k₁) made, then the likelihood factorizes and reduces to a product form

$\begin{array}{ll}{L}(A_2,A_1|\theta )=\left({\displaystyle \sum _{k_2=0}^n}\text{P}\left[A_2|k_2\right]\text{P}\left[k_2\right]\right)\left({\displaystyle \sum _{k_1=0}^n}\text{P}\left[A_1|k_1\right]\text{P}\left[k_1\right]\right)& (11)\end{array}$

that is much more computationally tractable in that only 2(n + 1) terms are required. This approach was taken by Costa et al. (2013) and combined with an additional approximation that neglected quantal synaptic components to focus on the mean effects of short-term plasticity. For the full probability density in which correlations are retained, it is not possible to factorize the likelihood into a scalar product in this way. However, we will show in the following sections that it is possible to use a Markovian property of this likelihood to factorize the calculation into a matrix product.

3.3. Joint Probability for Serial Release Events

The quantal component of the likelihood is most problematic; to illustrate the method of tractably calculating the full likelihood we will first consider the joint probability of paired release events P(k₂, k₁). The generalization to a train of many presynaptic action potentials is straightforward. Note that knowing the number of release events at a particular action potential does not specify the state of the system; however, knowing the number of occupied release sites before a spike does fully specify the state of system. This is the Markovian property that makes likelihood calculation possible. We call y_m the number of available vesicles present just before the mth action potential. Note that the expected value of y_m, $\text{E}\left[y_m\right]=n{\langle x\rangle }_m^{\ominus }$, where ${\langle x\rangle }_m^{\ominus }$ obeys Equation (2). Using this notation we can write the paired release probability in a more verbose form

$\begin{array}{ll}\text{P}(k_2,k_1)={\displaystyle \sum _{y_2=0}^n}{\displaystyle \sum _{y_1=0}^n}\text{P}(k_2,y_2,k_1,y_1).& (12)\end{array}$

It is now possible to factorize the probability on the right-hand-side of the above equation. First we use the product rule to expand as follows

$\begin{array}{ll}\text{P}(k_2,y_2,k_1,y_1)=\text{P}(k_2,y_2,k_1|y_1)\text{P}(y_1)& (13)\end{array}$

$\begin{array}{ll}=\text{P}(k_2,y_2|k_1,y_1)\text{P}(k_1|y_1)\text{P}(y_1)& (14)\end{array}$

$\begin{array}{l}=\text{P}(k_2|y_2,k_1,y_1)\text{P}(y_2|k_1,y_1)\\ \text{P}(k_1|y_1)\text{P}(y_1)& (15)\end{array}$

$\begin{array}{ll}=\text{P}(k_2|y_2)\text{P}(y_2|k_1,y_1)\text{P}(k_1|y_1)\text{P}(y_1)& (16)\end{array}$

where in the last step we have used the Markovian property of the occupancy variable. Note also that this is an iterative procedure, in which we can factorize the joint probability starting with the first action potential and then the second, that can be continued for joint probabilities that are comprised of an arbitrary number of spikes. For example, for the case of three action potentials it is only necessary to multiply the two-spike case by P(k₃|y₃)P(y₃|k₂, y₂) with the generalization to higher numbers of spike trains obvious. Inserting the final result in Equation (16) of this factorization into Equation (12) results in the following form for the two-spike case

$\begin{array}{ll}\text{P}(k_2,k_1)={\displaystyle \sum _{y_2=0}^n}{\displaystyle \sum _{y_1=0}^n}\left[\text{P}(k_2|y_2)\right]\left[\text{P}(y_2|k_1,y_1)\text{P}(k_1|y_1)\right]\left[\text{P}(y_1)\right]& (17)\end{array}$

where the square parentheses have been used to isolate components depending on k₂ or k₁ or neither. This form looks like an inner product and can be written in matrix-vector form (using bra-ket notation) as

$\begin{array}{ll}\text{P}(k_2,k_1)=\langle l_2\left|q_1\right|r_0\rangle & (18)\end{array}$

where 〈l₂| is a row vector dependent on k₂, q₁ is an (n + 1) by (n + 1) matrix dependent on k₁ and |r₀〉 is a column vector that comprises the initial conditions. Typically P(y₁) = δ_{y1, n} so that |r₀〉 has one non-zero entry to indicate that the synapse is initially fully stocked with vesicles. Note also that the case of three action potentials is straightforward

$\begin{array}{ll}\text{P}(k_3,k_2,k_1)=\langle l_3\left|q_2{q}_1\right|r_0\rangle & (19)\end{array}$

with obvious generalization to higher numbers of spikes. The joint release probability can therefore be reduced to matrix multiplication. The entries of the left row vector and matrices generally comprise two forms. The first form is simply the number of release events k_m chosen from the occupancy y_m, using the current probability of release u_m and is therefore binomial

$\begin{array}{ll}\text{P}(k_m|y_m)=\left(\begin{array}{c}{y}_m\\ k_m\end{array}\right)u_m^{k_m}{(1-u_m)}^{y_m-k_m}.\hfill & (20)\end{array}$

The second form gives the occupancy y_m+1 given k_m releases from an occupancy y_m at the previous action potential. This implies that there were n − y_m + k_m empty release sites just after the mth pulse. We require there to be n − y_m+1 empty sites just before the (m + 1)th pulse which means that y_m+1 − y_m + k_m sites were restocked. Let g_m be the restock probability of a single empty release site between time t_m and t_m+1

$\begin{array}{l}\text{P}(y_{m+1}|k_m,y_m)=\left(\begin{array}{c}n-y_m+k_m\\ y_{m+1}-y_m+k_m\end{array}\right)g_m^{y_{m+1}-y_m+k_m}\hfill \\ {\left(1-g_m\right)}^{n-y_{m+1}}\hfill & (21)\end{array}$

where this quantity depends on the time between spikes for the synaptic-dynamics model (and all other common synaptic models).

3.4. Joint Probability for Serial EPSP Amplitudes

We can now use the factorized form for the serial quantal release events to calculate the full likelihood, which is the joint probability density of seeing amplitudes A₁ and A₂ given the parameter set.

$\begin{array}{l}{L}(A_2,A_1|\theta )={\displaystyle \sum _{y_2}}{\displaystyle \sum _{y_1}}\left[{\displaystyle \sum _{k_2=0}^n}\text{P}\left[A_2|k_2\right]\text{P}\left[k_2|y_2\right]\right]\\ \left[{\displaystyle \sum _{k_1=0}^n}\text{P}\left[y_2|k_1,y_1\right]\text{P}\left[A_1|k_1\right]\text{P}\left[k_1|y_1\right]\right]\\ \left[\text{P}(y_1)\right].& (22)\end{array}$

The probabilities P[A₁|k₁] and P[A₂|k₂] for the observed amplitudes given that a certain number of vesicles were released are defined by Equation (8). The form of Equation (22) can again be interpreted as an inner product which can be written in bra-ket notation

$\begin{array}{ll}{L}(A_2,A_1|\theta )=\langle L_2\left|Q_1\right|R_0\rangle & (23)\end{array}$

where 〈L₂| is a row vector dependent on A₂, Q₁ is a matrix dependent on A₁ and |R₀〉 is a column vector with the initial configuration before the first action potential. This quantity is relatively straightforward to compute and, importantly, does not grow exponentially in computational complexity for higher numbers of action potentials. For example, for three spikes we have

$\begin{array}{ll}{L}(A_3,A_2,A_1|\theta )=\langle L_3\left|Q_2{Q}_1\right|R_0\rangle & (24)\end{array}$

with the generalization to higher numbers of presynaptic spikes straightforward.

3.5. Inferring Quantal Parameters from Synthetic Data

The methodology just described is first applied to synthetic data to test how well the correlated likelihood function can recover quantal and dynamic parameters (Figure 1). Here the synaptic-dynamics model is used to generate sweeps of synthetic amplitude trains. For this model, the eight parameters to infer are the release site number n, initial release probability p₀, facilitated release probability after an isolated spike p₁, depression timescale τ_D, facilitation timescale τ_F, mean quantal amplitude μ_a, standard deviation in quantal amplitude σ_a, and standard deviation of background noise σ_b.

FIGURE 1

Figure 1. Bayesian inference provides parameter distributions from five sweeps of synthetic data comprising 30 regular spikes at 30 Hz. Marginal posterior distributions (black), maximum a-posteriori estimates (orange crosses) and true parameter values (light blue dots) for the parameters of the synaptic model summarized in Table 1. Posteriors shown after 10⁶ Metropolis-Hastings samples. The true values were n = 7, τ_D = 0.25 s, τ_F = 0.2 s, p₀ = 0.6, p₁ = 0.8, μ_a = 0.25 mV, σ_a = 0.1 mV and σ_b = 0.05 mV.

Figure 1 shows marginal posterior distributions of these eight parameters given five simulated sweeps, each of 30 regular spikes at 30Hz. The posterior distributions reflect the true parameters well for all synaptic parameters with the exception of the facilitation timescale τ_F and quantal amplitude standard deviation σ_a. These parameters have been observed to be hard to estimate in previous studies, with Costa et al. (2013) finding broad distributions for τ_F, and Bhumbra and Beato (2013) and Barri et al. (2016) noting similar uncertainties in their estimates of quantal variability. The correlated Bayesian method does not qualitatively change these results, but makes the best use of available data to accurately estimate the uncertainty. The posterior distributions narrow with more data, but it is also possible to change experimental protocols to improve estimates. Costa et al. (2013) note that when the stimulation process is Poisson, rather than periodic, estimates of the time constants τ_D and τ_F using their method are improved due to the broader range of interspike intervals. This is equally true of the correlated Bayesian method. Estimates of σ_a could be improved by a very high stimulation rate that typically causes either 0 or 1 vesicles to release with each spike. Note that with typical delays between sweeps of 15 s, collecting this dataset required just over a minute of experimental time, giving a relatively sparse dataset that nevertheless still allows good estimates of the underlying synaptic parameters.

A major advantage of the Bayesian method over a maximum likelihood approach is that it can recover the full distribution of parameters. This allows determination of the covariances between different parameters. Figure 2 plots the joint posterior distributions of certain pairs of parameters (in total there are 28 possible pairs for the synaptic-dynamics model considered here). Figure 2A shows the relationship between release site number n, depression timescale τ_D, initial release probability p₀, and mean quantal amplitude μ_a. The inverse relationship between estimates of n and μ_a can be anticipated beause the mean EPSP size will always depend on the product of these two quantities. Note in particular that the relationship between release probability and both n and τ_D has a characteristic curved shape that is not apparent from looking at the individual marginal distributions. This is even more apparent (Figure 2B) for larger values of n that can be seen in some central synapses (Loebel et al., 2009, 2013).

FIGURE 2

Figure 2. Joint parameter estimates for the synaptic-dynamics model. (A) Pairwise and individual posterior marginals for release-site number n, depression timescale τ_D, initial release probability p₀, and mean quantal amplitude μ_a. True parameter values and data are the same as Figure 1. Colorbars for the values of the posterior distributions are not shown; the relative differences in value show the shape and sharpness of the pairwise posteriors for each pair of parameters. (B) Pairwise posterior marginal for release site number n and initial release probability p₀ for a case where the true values were n = 35 and p₀ = 0.50 showing a strong anticorrelation. All posteriors shown after 10⁶ Metropolis-Hastings samples.

3.6. Experiment: Changing Synaptic Dynamics under Adenosine Application

The neuromodulator adenosine is implicated (Kerr et al., 2013) in the developmental shift from dominant depression at juvenile synapses to weak facilitation at mature synapses (Reyes and Sakmann, 1999). Adenosine acts via A1 receptors to ultimately reduce the probability of vesicle release (Dunwiddie and Fredholm, 1989). Measurement of synaptic dynamics under control conditions and then during bath-application of adenosine therefore provides a convenient experimental protocol to test the inference method. For the control case an initially depressing juvenile connection was stimulated 40 times with nine periodic presynaptic spikes at 40 and 20 Hz (see Figure 3A) followed by a recovery spike, with the postsynaptic response recorded. Adenosine (100 μM) was then bath-applied to the slice (see Methods) and the stimulus protocol repeated.

FIGURE 3

Figure 3. Bayesian inference captures the shift in synaptic dynamics under application of adenosine. (A) Individual postsynaptic voltage traces under control (top) and adenosine (bottom) conditions. (B) Mean EPSP size for each spike in the stimulation protocol under control (blue) and adenosine (red) conditions. Bars show standard error. (C) Marginal posterior distributions for the parameters of the synaptic model in the control (blue) and adenosine (red) conditions. (D) Pairwise posterior marginals for number of active release sites n and initial release probability p₀ before (left) and after (right) application of adenosine. Posteriors shown after 5 × 10⁶ Metropolis-Hastings samples.

Figure 3A plots individual postsynaptic voltage traces before and after the application of adenosine; Figure 3B shows the change in average EPSP size. The marginal maximum-likelihood estimates for the depression timescale τ_D and mean quantal amplitude μ_a are similar between the control and adenosine datasets (Figure 3C). However, the suppressive effect of adenosine on synaptic transmission is clearly visible in the effective number of release sites n and the initial release probability p₀ that drives the shift from predominantly depressing to weakly facilitating synapses. It is also possible to examine the changes in covariance between pairs of parameters inferred from the experimental data (Figure 3D). Considering active release sites n and initial release probability p₀ together makes particularly apparent the shift in synaptic transmission.

3.7. Comparison with Methods that Neglect Serial Correlations

Previous Bayesian inference methods have demonstrated that an uncorrelated likelihood function can accurately infer the quantal (Bhumbra and Beato, 2013) and mean dynamic (Costa et al., 2013) parameters of a synapse. It can therefore be asked under what conditions does the exact likelihood function, which accounts for correlations, provide an improvement over existing methods. Synapses with low numbers of release sites n, high release probabilities u, or long depression timescales τ_D have the strongest correlations between EPSPs. High release probabilities u can arise either at strongly depressing synapses, with a high value of p₀, or facilitating synapses where the stimulation protocol causes large values of u(t) to arise. In addition to these, at least partly, physiological factors, the correlated likelihood function is superior in conditions of sparse data. When only a few PSPs are available per sweep or, more importantly, only a few sweeps are available correlations within a spike train are relatively more important. To quantify this, we compared the full likelihood function described above with an approximated likelihood calculated by ignoring correlations (calculated using forms like Equation 11). The approximate likelihood did not account for the observed previous PSP amplitudes within a sweep, only their distribution of probabilities given by the model parameters and previous spike times. As expected, the uncorrelated likelihood function gave broader posterior distributions (Figure 4A) with this effect diminishing as more data is added, either in the form of more EPSPs per sweep or more independent sweeps (Figure 4B). Overall, the exact likelihood function that accounts for correlations provides superior inference on synaptic parameters. It is possible to obtain accurate constraints on synaptic parameters with only a few sweeps, meaning that experiments could capture a snapshot of synaptic properties in a small time window during protocols that change synaptic properties on timescales of tens of seconds rather than tens of minutes.

FIGURE 4

Figure 4. Comparison of likelihood functions that do or do not account for serial correlations in synaptic amplitudes. (A) Posterior distributions for release site number n computed by correlated (solid) and uncorrelated (dashed) likelihood functions for three different values of n (n = 8 dark black; n = 15 middle blue; and n = 35 light blue) for a single sweep of 5 spikes regularly distributed at 30Hz. (B) 95% confidence intervals for correlated (solid) and uncorrelated (dashed) likelihood functions as a function of the number of sweeps for different numbers of spikes per train. Spikes occur at 30Hz, the true value of n is 35, and averages are taken over 10 realizations. Other parameters are the same as for Figure 1 (light-blue dots).

4. Discussion

We have presented a method for exactly and efficiently calculating the probability of a given train of PSP amplitudes for dynamical synapses with the utility and robustness of the method demonstrated on synthetic and experimental data. This method, presented earlier in Bird (2016) is equivalent to that simultaneously and independently discovered by Barri et al. (2016) in their expectation-maximization approach, and represents a combination and extension of the recent work of Bhumbra and Beato (2013) on the exact likehood of isolated events and Costa et al. (2013) on the approximated likelihood of serial events. By considering quantal and dynamic properties together, the method described accounts for information that is necessarily neglected when each component is examined in isolation. The advance renders the calculation of the likelihood required for Bayesian inference practical for a variety of models of short-term synaptic plasticity. Moreover, unlike approaches that have relied on mean-variance analysis, it is applicable to single-sweep experiments and so is suitable for in-vivo scenarios where presynaptic firing is uncontrolled, but can be monitored.

The likelihood calculation that makes this inference possible is flexible and can be extended to a number of common synaptic models, allowing for examination of augmented recovery (Wang and Kaczmarek, 1998; Hosoi et al., 2007), release-independent depression with frequency-dependent recovery (Fuhrmann et al., 2004), and receptor desensitization (Jones and Westbrook, 1996; Otis et al., 1996). Four such models are described in Appendix A with associated computer code in the MATLAB and JULIA environments to be found in the Supplementary Material. Another natural and straightforward extension of the methodology presented here is to not assume that all sites are initially occupied but have the initial state of the system as a parameter to be inferred. This scenario is relevant for in-vivo experiments where there is no natural break in the presynaptic activity: in this case the release site occupancy and state of the dynamic release probability would be unknown.

Author Contributions

AB and MR wrote the paper, AB and MR derived the equations, AB made the figures, AB wrote the MATLAB code and MR wrote the JULIA code, MW supervised experimental data collection.

Funding

We acknowledge funding through a Warwick Systems Biology Doctoral Training Centre fellowship to AB, funded by the UK BBSRC funding agency (BBSRC Grant No. BB/G530233/1), and funding to MR under BBSRC Grant No. BB/J015369/1.

Conflict of Interest Statement

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.

Acknowledgments

Data under adenosine was gathered by Dr. Michael Kerr. Thanks to Dr. Adam Newton (Warwick) and Dr. Peter Jedlicka (Frankfurt) for helpful comments on the manuscript. MR thanks Dr. Christophe Ladroue for useful discussions.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/article/10.3389/fncom.2016.00116/full#supplementary-material

Supplementary Material comprises MATLAB and JULIA code to run Bayesian inference using the depression (DEP), depression and facilitation (FAD), release-independent depression (RID) and release-independent depression with frequency-dependent recovery (FDR) synaptic models described, in the Methods and Appendix A.

References

Abbott, L. F. (1997). Synaptic depression and cortical gain control. Science 275, 221–224. doi: 10.1126/science.275.5297.221

PubMed Abstract | CrossRef Full Text | Google Scholar

Abbott, L. F., and Regehr, W. G (2004). Synaptic computation. Nature 431, 796–803. doi: 10.1038/nature03010

PubMed Abstract | CrossRef Full Text | Google Scholar

Barri, A., Wang, Y., Hansel, D., and Mongillo, G. (2016). Quantifying repetitive transmission at chemical synapses: a generative-model approach. eNeuro. 3:ENEURO.0113-15.2016. doi: 10.1523/ENEURO.0113-15.2016

PubMed Abstract | CrossRef Full Text | Google Scholar

Bekkers J. (1994). Quantal analysis of synaptic transmission in the central nervous system. Curr. Opin. Neurobiol. 4, 360–365. doi: 10.1016/0959-4388(94)90097-3

PubMed Abstract | CrossRef Full Text

Bennett, M. R., and Florin, T. (1974). Statistical analysis of the release of acetylcholine at newly formed synapses in striated muscle. J. Physiol. 238, 93–107. doi: 10.1113/jphysiol.1974.sp010512

PubMed Abstract | CrossRef Full Text | Google Scholar

Bennett, M. R., and Kearns, J. L. (2000). Statistics of transmitter release at nerve terminals. Prog. Neurobiol. 60, 545–606. doi: 10.1016/S0301-0082(99)00040-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhumbra, G. S., and Beato, M. (2013). Reliable evaluation of the quantal determinants of synaptic efficacy using Bayesian analysis. J. Neurophysiol. 109, 603–620. doi: 10.1152/jn.00528.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Bird, A. D. (2016). Temporal and Spatial Factors Affecting Synaptic Transmission in Cortex. Ph.D. Thesis, University of Warwick.

Birò, A., Holderith, N., and Nusser, Z. (2005). Quantal size is independent of the release probability at hippocampal excitatory synapses. J. Neurosci. 25, 223–232. doi: 10.1523/JNEUROSCI.3688-04.2005

PubMed Abstract | CrossRef Full Text

Boyd, I. A., and Martin, A. R. (1956). The end-plate potential in mammalian muscle. J. Physiol. 132, 74–91. doi: 10.1113/jphysiol.1956.sp005503

PubMed Abstract | CrossRef Full Text | Google Scholar

Brémaud, A., West, D., and Thomson, A. (2007). Binomial parameters differ across neocortical layers and with different classes of connections in adult rat and cat neocortex. Proc. Natl. Acad. Sci. U.S.A. 104, 14134–14139. doi: 10.1073/pnas.0705661104

PubMed Abstract | CrossRef Full Text | Google Scholar

Clements, J. D. (2003). Variance-mean analysis: a simple and reliable approach for investigating synaptic transmission and modulation. J. Neurosci. Meth. 130, 115–125. doi: 10.1016/j.jneumeth.2003.09.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Costa, R. P., Sjöström, P. J., and van Rossum, M. C. W. (2013). Probabilistic inference of short-term synaptic plasticity in neocortical microcircuits. Front. Comput. Neurosci. 7:75. doi: 10.3389/fncom.2013.00075

PubMed Abstract | CrossRef Full Text | Google Scholar

del Castillo, J., and Katz, B. (1954). Quantal components of the end-plate potential. J. Physiol. 124, 560–573. doi: 10.1113/jphysiol.1954.sp005129

PubMed Abstract | CrossRef Full Text | Google Scholar

Dudel, J., and Kuffler, S. (1961). Mechanism of facilitation at the crayfish neuromuscular junction. J. Physiol. 155, 530–542. doi: 10.1113/jphysiol.1961.sp006645

PubMed Abstract | CrossRef Full Text | Google Scholar

Dunwiddie, T. V., and Fredholm, B. B. (1989). Adenosine A1 receptors inhibit adenylate cyclase activity and neurotransmitter release and hyperpolarise pyramidal neurons in rat hippocampus. J. Pharmacol. Exp. Ther. 249, 31–37.

PubMed Abstract | Google Scholar

Eccles, J. C., Katz, B., and Kuffler, S. W. (1941). Nature of the endplate potential in curarized muscle. J. Neurophysiol. 4, 362–387.

Google Scholar

Fatt, P., and Katz, B. (1952). Spontaneous subthreshold activity at motor nerve endings. J. Physiol. 117, 109–128.

PubMed Abstract | Google Scholar

Foster, K., and Regehr, W. (2004). Variance-mean analysis in the presence of a rapid antagonist indicates vesicle depletion underlies depression at the climbing fiber synapse. Neuron 43, 119–131. doi: 10.1016/j.neuron.2004.06.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Franks, K., Stevens, C., and Sejnowski, T. (2003). Independent sources of quantal variability at single glutamatergic synapses. J. Neurosci. 23, 3186–3195.

PubMed Abstract | Google Scholar

Fuhrmann, G., Cowan, A., Segev, I., Tsodyks, M., and Stricker, C. (2004). Multiple mechanisms govern the dynamics of depression at neocortical synapses of young rats. J. Physiol. 557, 415–438. doi: 10.1113/jphysiol.2003.058107

PubMed Abstract | CrossRef Full Text | Google Scholar

Fuhrmann, G., Segev, I., Markram, H., and Tsodyks, M. (2002). Coding of temporal information by activity-dependent synapses. J. Neurophysiol. 87, 140–148. doi: 10.1152/jn.00258.2001

PubMed Abstract | CrossRef Full Text | Google Scholar

Furukawa, T., Kuno, M., and Matsuura, S. (1982). Quantal analysis of a decremental response at hair cell-afferent fibre synapses in the goldfish sacculus. J. Physiol. 322, 181–195. doi: 10.1113/jphysiol.1982.sp014031

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanse, E., and Gustafsson, B. (2001). Quantal variability at glutamatergic synapses in area CA1 of the rat neonatal hippocampus. J. Physiol. 531, 467–480. doi: 10.1111/j.1469-7793.2001.0467i.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hardingham, N. R., Read, J. C. A., Trevelyan, A. J., Nelson, J. C., Jack, J. J. B., and Bannister, N. J. (2010). Quantal analysis reveals a functional correlation between presynaptic and postsynaptic efficacy in excitiatory connections from rat neocortex. J. Neurosci. 30, 1441–1451. doi: 10.1523/JNEUROSCI.3244-09.2010

PubMed Abstract | CrossRef Full Text

Hastings, W. K. (1970). Monte Carlo sampling methods using Markov chains and their applications. Biometrika 57, 97–109. doi: 10.1093/biomet/57.1.97

CrossRef Full Text

Hosoi, N., Sakaba, T., and Neher, E. (2007). Quantitative analysis of calcium-dependent vesicle recruitment and its functional role at the Calyx of Held synapse. J. Neurosci. 27, 14286–14298. doi: 10.1523/JNEUROSCI.4122-07.2007

PubMed Abstract | CrossRef Full Text

Jones, M. V., and Westbrook, G. L. (1996). The impact of receptor desensitization on fast synaptic transmission. Trends Neurosci. 19, 96–101. doi: 10.1016/S0166-2236(96)80037-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Kerr, M. I., Wall, M. J., and Richardson, M. J. E. (2013). Adenosine A₁ receptor activation mediates the developmental shift at layer 5 pyramidal cell synapses and is a determinant of mature synaptic strength. J. Physiol. 591(Pt 13), 3371–3380. doi: 10.1113/jphysiol.2012.244392

PubMed Abstract | CrossRef Full Text | Google Scholar

Korn, H., Triller, A., Mallet, A., and Faber, D. S. (1981). Fluctuating responses at a central synapse: n of binomial fit predicts number of stained presynaptic boutons. Science 213, 898–901. doi: 10.1126/science.6266015

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuno, M. (1964). Quantal components of excitatory synaptic potentials in spinal motoneurones. J. Physiol. 175, 81–99. doi: 10.1113/jphysiol.1964.sp007504

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuno, M. (1971). Quantum aspects of central and ganglionic synaptic transmission in vertebrates. Physiol. Rev. 51, 647–678.

PubMed Abstract | Google Scholar

Kuno, M., and Weakly, J. N. (1972). Quantal components of the inhibitory synaptic potential in spinal mononeurones of the cat. J. Physiol. 224, 287–303. doi: 10.1113/jphysiol.1972.sp009895

PubMed Abstract | CrossRef Full Text | Google Scholar

Lefort, S., Tomm, C., Sarria, J. C. F., and Petersen, C. C. H. (2009). The excitatory neuronal network of the C2 barrel column in mouse primary somatosensory cortex. Neuron 61, 301–316. doi: 10.1016/j.neuron.2008.12.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Liley, A. W. (1956). The quantal components of the mammalian end-plate potential. J. Physiol. 133, 571–587. doi: 10.1113/jphysiol.1956.sp005610

PubMed Abstract | CrossRef Full Text | Google Scholar

Loebel, A., Le Bè, J. V., Richardson, M. J., Markram, H., and Herz, A. V. (2013). Matched pre- and post-synaptic changes underlie synaptic plasticity over long time scales. J. Neurosci. 33, 6257–6266. doi: 10.1523/JNEUROSCI.3740-12.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Loebel, A., Silberberg, G., Helbig, D., Markram, H., Tsodyks, M., and Richardson, M. J. E. (2009). Multiquantal release underlies the distribution of synaptic efficacies in the neocortex. Front. Comput. Neurosci. 3:27. doi: 10.3389/neuro.10.027.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

McGuinness, L., Taylor, C., Taylor, R. D. T., Yau, C., Langenhan, T., Hart, M. L., et al. (2010). Presynaptic NMDARs in the hippocampus facilitate transmitter release at theta frequency. Neuron 68, 1109–1127. doi: 10.1016/j.neuron.2010.11.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., and Teller, E. (1953). Equation of state calculations by fast computing machines. J. Chem. Phys. 21, 1087–1092. doi: 10.1063/1.1699114

CrossRef Full Text | Google Scholar

Mongillo, G., Barak, O., and Tsodyks, M. (2008). Synaptic theory of working memory. Science 319, 1543–1546. doi: 10.1126/science.1150769

PubMed Abstract | CrossRef Full Text | Google Scholar

Otis, T., Zhang, S., and Trussel, L. O. (1996). Direct measurement of AMPA receptor desensitization induced by glutamatergic synaptic transmission. J. Neurosci. 16, 7496–7504.

PubMed Abstract

Reyes, A., and Sakmann, B. (1999). Developmental switch in the short-term modification of unitary EPSPs evoked in layer 2/3 and layer 5 pyramidal neurons of rat neocortex. J. Neurosci. 19, 3827–3835.

PubMed Abstract | Google Scholar

Richardson, M. J. E., and Silberberg, G. (2008). Measurement and analysis of postsynaptic potentials using a novel voltage-deconvolution method. J. Neurophysiol. 99, 1020–1031. doi: 10.1152/jn.00942.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, J. (1976). Estimation of parameters for a model of transmitter release at synapses. Biometrics 32, 61–68. doi: 10.2307/2529338

PubMed Abstract | CrossRef Full Text | Google Scholar

Silver, R. (2003). Estimation of nonuniform quantal parameters with multiple-probability fluctuation analysis: theory, application and limitations. J. Neurolsci. Meth. 130, 127–141. doi: 10.1016/j.jneumeth.2003.09.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Silver, R. A., Momiyama, A., and Cull-Candy, S. (1998). Locus of frequency- dependent depression identified with multiple-probability fluctuation analysis at rat climbing fibre-Purkinje cell synapses. J. Physiol. 510, 881–902. doi: 10.1111/j.1469-7793.1998.881bj.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, S., Sjöström, P. J., Reigl, M., Nelson, S., and Chklovskii, D. B. (2005). Highly nonrandom features of synaptic connectivity in local cortical circuits. PLoS Biol. 3:68. doi: 10.1371/journal.pbio.0030068

PubMed Abstract | CrossRef Full Text | Google Scholar

Stein, R. B. (1965). A theoretical analysis of neuronal variability. Biophys. J. 5, 173–194. doi: 10.1016/S0006-3495(65)86709-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Thomson, A. M., Deuchars, J., and West, D. C. (1993). Large, deep layer pyramid-pyramid single axon EPSPs in slices of rat motor cortex display paired pulse and frequency-dependent depression, mediated presynaptically and self-facilitation, mediated postsynaptically. J. Neurophysiol. 70, 2354–2369.

PubMed Abstract | Google Scholar

Tsodyks, M. V., and Markram, H. (1997). The neural code between neocortical pyramidal neurons depends on neurotransmitter release probability. Proc. Natl. Acad. Sci. U.S.A. 94, 719–723. doi: 10.1073/pnas.94.2.719

PubMed Abstract | CrossRef Full Text | Google Scholar

Tsodyks, M. V., Pawelzik, K., and Markram, H. (1998). Neural networks with dynamic synapses. Neural Comp. 10, 821–835. doi: 10.1162/089976698300017502

PubMed Abstract | CrossRef Full Text | Google Scholar

Turner, D., and West, M. (1993). Bayesian analysis of mixtures applied to post-synaptic potential fluctuations. J. Neurosci. Meth. 47, 1–21. doi: 10.1016/0165-0270(93)90017-L

PubMed Abstract | CrossRef Full Text | Google Scholar

Varela, J. A., Sen, K., Gibson, J., Fost, J., Abbott, L. F., and Nelson, S. B. (1997). A quantitative description of short-term plasticity at excitatory synapses in layer 2/3 of rat primary visual cortex. J. Neurosci. 17, 7926–7940.

PubMed Abstract | Google Scholar

Wang, L., and Kaczmarek, L. (1998). High-frequency firing helps replenish the readily releasable pool of synaptic vesicles. Nature 394, 384–388. doi: 10.1038/28645

PubMed Abstract | CrossRef Full Text | Google Scholar

Zucker, R. S., and Regehr, W. G. (2002). Short-term synaptic plasticity. Annu. Rev. Physiol. 64, 355–405. doi: 10.1146/annurev.physiol.64.092501.114547

PubMed Abstract | CrossRef Full Text | Google Scholar

Appendix A

Extension to Other Synaptic Models

The likelihood calculation that was illustrated in the main text for a model with depression and facilitation can be straightforwardly adapted to other commonly used models of synaptic dynamics. These comprise models in which the restock probability g_m between presynaptic action potentials m and m + 1 and probability of release at the arrival of the mth action potential u_m depends only on the pattern of presynaptic activity. As part of the Supplementary Material we provide computer code for four such models, which are now described below with the synaptic parameters and dynamic variables tabulated in Table A1.

TABLE A1

Table A1. Extended table of inferred parameters (top) and dynamic variables (bottom) used in the synaptic models discussed in Appendix A.

(i) Depression only - DEP model

This is perhaps the simplest model of short-term synaptic plasticity and features only vesicle depletion and restock (Tsodyks and Markram, 1997; Fuhrmann et al., 2002). The occupation of a single release site is governed by the stochastic differential Equation (1). The mean-occupancy recursion relations for the model are given by Equations (2) and (3) with a constant release probability u_m = p₀. The Poissonian restock of empty release sites occurs at a constant rate 1/τ_D and so in this case the restock probability g_m is given by Equation (4).

(ii) Depression and facilitation - DAF model

This is the model described in the main text (Varela et al., 1997; Tsodyks et al., 1998; Fuhrmann et al., 2002) and applies to facilitating synapses. The probability of restock is defined by Equation (4) and the probability of release u_m is given by recursion Equation (6).

(iii) Release-independent depression - RID model

This model was introduced (Fuhrmann et al., 2004) for synapses that do not display facilitation and considers a different form of depression which is uncorrelated with the preceding EPSP amplitudes. Release-independent depression is a reduction in release probability u_m caused by spiking activity which decays on a timescale τ_I0 (it can be thought of as a kind of anti-facilitation). The release probability immediately after an isolated pulse is again called p₁ but in contrast to facilitation p₁ < p₀. In this formalism the release probability u(t) obeys

$\begin{array}{ll}\frac{du}{dt}=\frac{p_0-u}{{\tau }_{I_0}}-\frac{u}{p_0}(p_0-p_1){\displaystyle \sum _{t_m}}\delta (t-t_m)& (\text{A}1)\end{array}$

where t_m are the times of the presynaptic action-potentials. The values of u(t) just after the mth and before the (m + 1)th action potentials ($u_m^{\oplus }$ and u_m+1 respectively) are defined by the following recursion relations

$\begin{array}{l}{u}_m^{\oplus }=u_m-u_m\left(\frac{p_0-p_1}{p_0}\right)\mathrm{\text{and}}\\ u_{m+1}=p_0+(u_m^{\oplus }-p_0)e^{-T_m/{\tau }_{I_0}}.& (\text{A}2)\end{array}$

The restock probability g_m is given by Equation (4) and is common to the previous two models.

(iv) Release-independent depression with frequency dependent recovery - FDR model

The recovery from release-independent depression is often seen to be frequency dependent (Fuhrmann et al., 2004). To account for this the timescale τ_I(t) is now a dynamic variable with initial value τ_I0, magnitude after an isolated spike of τ_I1 and decay timescale ς_I. The relevant equations are now

$\begin{array}{ll}\frac{du}{dt}=\frac{p_0-u}{{\tau }_I}-\frac{u}{p_0}(p_0-p_1){\displaystyle \sum _{t_m}}\delta (t-t_m)& (\text{A}3)\end{array}$

$\begin{array}{ll}\frac{d{\tau }_I}{dt}=\frac{{\tau }_{I_0}-{\tau }_I}{{\varsigma }_I}-\frac{{\tau }_I}{{\tau }_{I_0}}({\tau }_{I_0}-{\tau }_{I_1}){\displaystyle \sum _{t_m}}\delta (t-t_m)& (\text{A}4)\end{array}$

The dynamic RID recovery timescale τ_I(t) obeys Equation (A4). The values of τ_I(t) just after the mth and before the (m + 1)th action potentials (${\tau }_{Im}^{\oplus }$ and ${\tau }_{Im+1}^{\ominus }$ respectively) are defined by the following recursion relations

$\begin{array}{l}{\tau }_{Im}^{\oplus }={\tau }_{Im}^{\ominus }-{\tau }_{Im}^{\ominus }\left(\frac{{\tau }_{I_0}-{\tau }_{I_1}}{{\tau }_{I_0}}\right)\mathrm{\text{and}}\\ {\tau }_{Im+1}^{\ominus }={\tau }_{I_0}+({\tau }_{Im}^{\oplus }-{\tau }_{I_0})e^{-T_m/{\varsigma }_I}.& (\text{A}5)\end{array}$

The release probability u(t) obeys Equation (A3) and can also be defined recursively, with u_m+1 the release probability at the (m + 1)th spike at time t_m+1 given by

$\begin{array}{ll}{u}_{m+1}=p_0+\frac{(u_m^{\oplus }-p_0){({\tau }_{Im}^{\oplus })}^{{}^{\frac{{\varsigma }_I}{{\tau }_{I_0}}}}}{{\left[{\tau }_{Im}^{\oplus }-{\tau }_{I_0}\left(1-e^{\frac{T_m}{{\varsigma }_I}}\right)\right]}^{{}^{\frac{{\varsigma }_I}{{\tau }_{I_0}}}}}& (\text{A}6)\end{array}$

where $u_m^{\oplus }=u_m-u_m\left(\frac{p_0-p_1}{p_0}\right)$ is the release probability immediately following the mth spike.

Appendix B

The Likelihood Convolution Integral

The convolution integral for the amplitude distribution Equation (8) is computed a large number of times. There are two difficulties in doing this efficiently: (i) evaluating a gamma function for large shape parameters and (ii) finding reasonable bounds for the range of integration.

Evaluating Gamma Functions

When μ_a ≪ σ_a the argument of the gamma function in the denominator of Equation (8) can grow very large in order to normalize the distribution. To avoid issues with this, we note that Stirling's approximation allows evaluation of the gamma function with large arguments

$\begin{array}{l}\Gamma (n+1)\\ \approx \sqrt{2\pi n}{\left(\frac{n}{e}\right)}^n\left(1+\frac{1}{12n}+\frac{1}{288n^2}-\frac{139}{51840n^3}-\frac{571}{2488320n^4}\right)\end{array}$

and introduce a variable κ(n) such that

$\begin{array}{l}\kappa (n)=\frac{1}{\Gamma (n+1)}\sqrt{2\pi n}{\left(\frac{n}{e}\right)}^n\\ \approx {\left(1+\frac{1}{12n}+\frac{1}{288n^2}-\frac{139}{51840n^3}-\frac{571}{2488320n^4}\right)}^{-1}.\end{array}$

For small values of n the function κ(n) can be evaluated exactly, whereas for larger arguments the second, approximate form is used.

Bounding the Range of Integration

The second difficulty involves finding bounds for the range of integration in Equation (8). Introducing γ = β − 1 and y_* = γ/λ, and substituting κ(γ) for the gamma function, Equation (8) can be rewritten as

$\begin{array}{l}\text{P}\left[A|k\right]=\frac{\kappa (\gamma )}{2\pi {\sigma }_b}\frac{1}{\sqrt{y_{*}/\lambda }}{\displaystyle {\int }_0^{\infty }}dy\\ exp\left(-\lambda (y-y_{*})+\gamma log\left(y/y_{*}\right)-\frac{{(A-y)}^2}{2{\sigma }_b^2}\right)\end{array}$

which is in a form that is straightforward to integrate numerically, once reasonable integration bounds are identified. To this end we note that by completing the squares in the exponential and introducing z = y/σ_b the integrand (ignoring constant prefactors) can be written as

$\begin{array}{l}f(z)=z^{\gamma }{e}^{-{(z-2t)}^2/2}\end{array}$

where t = (A/σ_b − λσ_b)/2. The integrand has a single positive maximum $z_{*}=t+\sqrt{t^2+\gamma }$ (note that γ = β − 1 > 0 because the assumption is made that μ_a > σ_a, see Methods). Our method then is to perform the y integral over an equivalent range of z for which $f(z)>f(z_{*})e^{-x}$ where x is some sufficiently large power (typically in the range 10–15) with the boundaries of the z range efficiently found using a Newton-Raphson method.