# Bayesian Inference of Synaptic Quantal Parameters from Correlated Vesicle Release

^{1}Theoretical Neuroscience Group, Warwick Systems Biology Centre, University of Warwick, Coventry, UK^{2}Ernst Strüngmann Institute for Neuroscience, Max Planck Society, Frankfurt, Germany^{3}Frankfurt Institute for Advanced Studies, Frankfurt, Germany^{4}School of Life Sciences, University of Warwick, Coventry, UK^{5}Warwick Mathematics Institute, University of Warwick, Coventry, UK

Synaptic transmission is both history-dependent and stochastic, resulting in varying responses to presentations of the same presynaptic stimulus. This complicates attempts to infer synaptic parameters and has led to the proposal of a number of different strategies for their quantification. Recently Bayesian approaches have been applied to make more efficient use of the data collected in paired intracellular recordings. Methods have been developed that either provide a complete model of the distribution of amplitudes for isolated responses or approximate the amplitude distributions of a train of post-synaptic potentials, with correct short-term synaptic dynamics but neglecting correlations. In both cases the methods provided significantly improved inference of model parameters as compared to existing mean-variance fitting approaches. However, for synapses with high release probability, low vesicle number or relatively low restock rate and for data in which only one or few repeats of the same pattern are available, correlations between serial events can allow for the extraction of significantly more information from experiment: a more complete Bayesian approach would take this into account also. This has not been possible previously because of the technical difficulty in calculating the likelihood of amplitudes seen in correlated post-synaptic potential trains; however, recent theoretical advances have now rendered the likelihood calculation tractable for a broad class of synaptic dynamics models. Here we present a compact mathematical form for the likelihood in terms of a matrix product and demonstrate how marginals of the posterior provide information on covariance of parameter distributions. The associated computer code for Bayesian parameter inference for a variety of models of synaptic dynamics is provided in the Supplementary Material allowing for quantal and dynamical parameters to be readily inferred from experimental data sets.

## 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^{2+} 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^{2+} 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^{2+} 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^{2+} 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 *m*th 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

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

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 *m*th spike. Similarly, the probability of occupancy just before (*m* + 1)th AP can be related to the occupancy just after the *m*th AP as

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

where *T*_{m} = *t*_{m+1} − *t*_{m} is the time between the *m*th and (*m* + 1)th APs. Let *p*_{0} be the baseline value of the probability of release, and *p*_{1} 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*_{0} with timescale τ_{F}. The dynamics of *u*(*t*) therefore obeys

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

where the initial conditions are that *u*_{1} = *p*_{0} 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*_{1} = 2*p*_{0} (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

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

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 *m*th pulse then it is released with probability *u*_{m}. Between the *m*th 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^{2+} 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*_{0} of release for an isolated pulse; immediately after an isolated presynaptic action potential the release probability increases to *p*_{1}. The release probability *u* returns to its initial value *p*_{0} 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 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*_{1}, *t*_{2}, ⋯ , *t*_{M} and post-synaptic amplitudes *A*_{1}, *A*_{2} ⋯ *A*_{M}. The aim of the inference procedure is to calculate the probability densities of the parameters of the model θ = {*n, p*_{0}, *p*_{1}, τ_{D}⋯ } given the observed presynaptic action potential times {*t*_{1}, ⋯ , *t*_{M}} and postsynaptic amplitudes {*A*_{1}, ⋯*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:

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

where *k*_{m} is the number of vesicles released by the *m*th spike. Because of the nested sums there are (*n* + 1)^{2} 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*_{2}, *k*_{1}]; 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*_{2}, *k*_{1}) ≃ P(*k*_{2})*P*(*k*_{1}) made, then the likelihood factorizes and reduces to a product form

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*_{2}, *k*_{1}). 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 *m*th 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

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

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*_{3}|*y*_{3})P(*y*_{3}|*k*_{2}, *y*_{2}) 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

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

where 〈*l*_{2}| is a row vector dependent on *k*_{2}, *q*_{1} is an (*n* + 1) by (*n* + 1) matrix dependent on *k*_{1} and |*r*_{0}〉 is a column vector that comprises the initial conditions. Typically P(*y*_{1}) = δ_{y1, n} so that |*r*_{0}〉 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

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

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 *m*th 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}

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*_{1} and *A*_{2} given the parameter set.

The probabilities P[*A*_{1}|*k*_{1}] and P[*A*_{2}|*k*_{2}] 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

where 〈*L*_{2}| is a row vector dependent on *A*_{2}, *Q*_{1} is a matrix dependent on *A*_{1} and |*R*_{0}〉 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

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*_{0}, facilitated release probability after an isolated spike *p*_{1}, 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. 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^{6} Metropolis-Hastings samples. The true values were *n* = 7, τ_{D} = 0.25 s, τ_{F} = 0.2 s, *p*_{0} = 0.6, *p*_{1} = 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*_{0}, 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. 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*_{0}, 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*_{0} for a case where the true values were *n* = 35 and *p*_{0} = 0.50 showing a strong anticorrelation. All posteriors shown after 10^{6} 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 *A*1 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. 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*_{0} before (left) and after (right) application of adenosine. Posteriors shown after 5 × 10^{6} 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*_{0} 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*_{0} 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*_{0}, 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. 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

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

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

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

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

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

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

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

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

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

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

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

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

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

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.

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

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

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

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

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

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

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

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

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

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

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

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

Kerr, M. I., Wall, M. J., and Richardson, M. J. E. (2013). Adenosine A_{1} 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

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

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

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

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

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

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

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

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

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

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

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

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.

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.

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

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

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

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

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

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

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.

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

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

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

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.

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

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

## 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 *m*th 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. 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*_{0}. 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*_{1} but in contrast to facilitation *p*_{1} < *p*_{0}. In this formalism the release probability *u*(*t*) obeys

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

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

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

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

where ${u}_{m}^{\oplus}={u}_{m}-{u}_{m}\left(\frac{{p}_{0}-{p}_{1}}{{p}_{0}}\right)$ is the release probability immediately following the *m*th 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

and introduce a variable κ(*n*) such that

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

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

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.

Keywords: correlation, Bayesian, EPSP, synapse, quantal, stochastic, plasticity

Citation: Bird AD, Wall MJ and Richardson MJE (2016) Bayesian Inference of Synaptic Quantal Parameters from Correlated Vesicle Release. *Front. Comput. Neurosci*. 10:116. doi: 10.3389/fncom.2016.00116

Received: 26 August 2016; Accepted: 28 October 2016;

Published: 25 November 2016.

Edited by:

David Holcman, École Normale Supérieure, FranceReviewed by:

Gianluigi Mongillo, Paris Descartes University, FranceMark C. W. Van Rossum, University of Edinburgh, UK

Copyright © 2016 Bird, Wall and Richardson. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Alex D. Bird, a.d.bird@warwick.ac.uk