Parsing Out the Variability of Transmission at Central Synapses Using Optical Quantal Analysis

1. Introduction

Our understanding of the factors that contribute to the stochastic and variable process of synaptic transmission has improved steadily over the last few decades (Branco and Staras, 2009; Ribrault et al., 2011; Llera-Montero et al., 2019). It is now generally agreed that, at most glutamatergic synapses, quantal release does not saturate postsynaptic receptors (Liu et al., 1999, 2003; Umemiya et al., 1999; McAllister and Stevens, 2000; Nimchinsky et al., 2004) and that variability in trial-to-trial neurotransmission arises primarily from differences in the profile of glutamate released into the synaptic cleft (Ribrault et al., 2011). Several presynaptic mechanisms have been proposed to account for such amplitude fluctuations – uneven packaging of glutamate into synaptic vesicles, differences in release location within a synaptic terminal, diffusion process in the synaptic cleft and mode of exocytosis (Choi et al., 2003; Franks et al., 2003; Wu et al., 2007; Richards, 2009). As an additional factor, multiquantal release has been observed at many central synapses (Auger et al., 1998; Wadiche and Jahr, 2001; Oertner et al., 2002; Conti and Lisman, 2003; Christie and Jahr, 2006; Higley et al., 2009; Rudolph et al., 2015), where two or more vesicles are released quasi simultaneously at single synapses in response to the same electrical stimulus. Since each of these sources of variability impact the transmission of information differently, it is therefore important to parse out the relative proportion of different sources of variability at central synapses.

Several experimental methodologies have been developed to monitor transmission at single synapses (Malinow and Tsien, 1990; Dobrunz and Stevens, 1997, 1999; Oertner et al., 2002; Rudolph et al., 2015). Here we describe an optical-based technique and provide a number of validation and calibration experiments for the intensity-based optical glutamate sensor, iGluSnFR (Marvin et al., 2013), for optical quantal analysis at central synapses. We further provide a detailed theoretical and quantitative analysis for estimating fundamental features of quantal glutamate release. Leveraging experimental and statistical techniques, combined with a theoretically sound model, we present a formalism that is well poised to parse out the structure of variability and information content at central synapses.

2. Results

To study quantal features of glutamate release at central synapses, we turned to a genetically encoded intensity-based glutamate sensing florescent reporter (iGluSnFR). The versatility and usefulness of iGluSnFR – a diffuse, plasma membrane-bound optical reporter of glutamate release – has been demonstrated in both microscopic and macroscopic brain compartments (Borghuis et al., 2013; Marvin et al., 2013; Park et al., 2014; Parsons et al., 2016; Xie et al., 2016), although it has relatively seldom been used to study features of glutamate release at single spines (Soares et al., 2017; Jensen et al., 2019). To this end, we introduced iGluSnFR along with the morphological marker mCherry to CA1 neurons in hippocampal organotypic slices using biolistic transfection several days prior the experiments (Figure 1A; Soares et al., 2014, 2017). A detailed description of these procedures is available in Soares et al. (2014). This overall approach was favored since it allows for sparse transfection thereby allowing us to resolve optical signals from single spines with high contrast. Transfected neurons were imaged by two-photon microscopy using an excitation wavelength of 950 nm (Figure 1B) which we found to allow detection of both the iGluSnFR and mCherry fluorescent signal simultaneously. Dendritic spines in the apical dendritic arbor of transfected CA1 neurons were targeted for optical quantal analysis experiments. These contacts are likely the postsynaptic targets of Schaffer's collateral axons.

2.1. iGluSnFR-Mediated Monitoring of Endogenous Glutamate Release

Pyramidal neurons were identified by their localization in the slice and morphology. Namely imaging targeted to the CA1 region and we sought the clear presence of basal and apical spinous dendritic arborisation. The morphological identification was typically carried out by solely monitoring mCherry fluorescence. However, the baseline iGluSnFR fluorescence was typically fairly high, homogenously distributed across neuronal compartments and spines were readily observable, thereby readily allowing for broad cell-type identification. A typical experiment began by randomly surveying the apical arbor of an iGluSnFR-expressing cell for dendritic spines that exhibit a time-locked fluorescent responses to electrical stimuli delivered via a glass pipette positioned in stratum radiatum. In a few experiments, Alexa 594 was included in the internal solution of the stimulating electrode for direct visualization (Figure 1B), however, in the majority of experiments this dye was omitted and the stimulating electrode was maneuvered in the slice under visual guidance solely using differential interference contrast microscopy.

The optical detection of synaptic events that are eminently short-lived, spatially distributed and scarce is inherently challenging and deserves attention. In principle, imaging in frame scanning mode would be ideal to monitor synaptic fluorescent events from large dendritic regions, but it is hindered by limited signal to noise ratio and temporal resolution (Figure 1C). We thus carried out line scan experiments wherein multiple neighboring spines were monitored simultaneously (Figure 1D). This approach offered the ability to survey multiple spines at once with a scanning frequency (>500 Hz) sufficient to visually identify rapid glutamate transients. To circumvent the relative paucity of synaptic events due to the probabilistic nature of release, paired-pulse electrical stimulation (50-100 ms inter-stimulus interval) were delivered to increase the likelihood of release during this initial probing phase. Lastly, realistic range of stimulus intensities was determined by parallel and historical whole-cell electrophysiological recordings by the same experimenter. Once a responsive spine was identified, a short line scan was redrawn through the spine and its parent dendritic compartment to capture the spatial profile of glutamate release. The electrical stimulation was then gradually reduced to the minimal intensity that still evoked time-locked responsiveness. This last step was taken in order to reduce the potential of signal contamination by glutamate spillover from neighboring synapses. The identified spines routinely stayed responsive to electrical stimulation for long durations (> 1 hour), opening the door to the repetitive low frequency sampling methodology required for building a dataset sufficient for optical quantal analysis.

2.2. Extraction of Regions of Interest

Spatial discrimination of iGluSnFR signals emanating from either spine or dendritic compartments was achieved by analyzing the intensity profile across the line scan, which was drawn orthogonal to the parent dendrite. The trough between spine and dendrite peaks was used to split the signal of the line scan into the two compartments (Figure 1E) to isolate spine- and dendrite-specific iGluSnFR transients (Figure 1D, right). Larger amplitude iGluSnFR transients were generally observed in the spine compartment, indicating that the density of glutamate release was mostly concentrated at the spine. When present, the dendritic fluorescence transients were of smaller amplitudes and co-varied with that recorded from the spine compartment, suggesting that dendritic signals were likely the result of spillover from the parent spine rather than from release from a distinct, neighboring synapse (Figure 1F). As such, we used only the spine compartment signal for all subsequent analyses. Finally, and consistent with the probabilistic release of glutamate vesicles at these synapses, release failures were readily observed (Figure 1G). These results demonstrate that iGluSnFR is a useful optical reporter for single-spine quantal analysis.

2.3. Glutamate vs. Post-synaptic Calcium Sensors for Opto-Quantal Analysis

The difficulty in unambiguously and routinely study neurotransmitter release from a single synapse due to the lack of spatial resolution afforded by electrophysiological recordings has been a longstanding shortcoming. By providing spatial information, optically-based approaches for quantal analysis offers promise of a solution to this problem, yet are limited by temporal resolution generally poorer than that afforded by cellular electrophysiology. By using two-photon uncaging of MNI-glutamate to precisely control the amount and timing of glutamate released onto single spines, we next sought to examine the kinetic performances of iGluSnFR by a side-by-side comparison with other commonly used reporters of glutamate release for quantal analysis. Specifically, we sought to compare with electrophysiological monitoring of synaptic AMPAR activation and optical recordings of quantal analysis using NMDAR-mediated calcium influx by calcium indicators. Since optical recordings of calcium influx using the GCaMP family of genetically encoded calcium-indicators are becoming increasingly popular, we turned our attention to GCaMP6f, a fast variant of the GCaMP family.

We obtained whole-cell recordings from CA1 neurons transfected with either iGluSnFR of GCaMP6f (Figure 2A) and voltage-clamped the cell at -70 mV. While continuously imaging the spine of interest (at ≈ 715 Hz), a second laser line tuned to 720 nm delivered a 1 ms light pulse to the tip of the spine in the presence of MNI-Glutamate (2.5mM), to induce uncaging-evoked optical transients recorded simultaneously with EPSCs (Figures 2B,C). In response to repetitive presentation of nominally constant concentration of glutamate by 2P uncaging at single synapses (Figure 2D), we compared the performance of 3 distinct reporters of glutamate transients at single synapses: (i) iGluSnFR transients; (ii) GCaMP6f transients (i.e., NMDAR-dependent calcium influx) and, (iii) AMPAR activation (uncaging-evoked EPSCs; uEPSCs). We found that the trial-to-trial variability of the iGluSnFR responses was remarkably low, even lower than that of uEPSCs (Figure 2E). In keeping with the more complex and convolved nature of NMDAR- and calcium-mediated optical detection of glutamate release, GCaMP6f-mediated signal displayed the largest variability of the 3 approaches (Figure 2E). iGluSnFR transients also displayed much faster decay kinetics (Figure 2F) and rise time (Figure 2G) compared to GCaMP6f, and were remarkably close but still slower than the kinetics of uEPSCs (p < 0.001 for both comparisons, unpaired student's t-test). The kinetic properties of iGluSnFR in response to glutamate uncaging therefore favorably compares to those of the calcium-sensitive organic dyes Alexa 4FF (Lee et al., 2016) and Oregon Green BAPTA-1 (unpublished observations) that are significantly slower. The fast kinetics of iGluSnFR enable discrimination of successive stimulus peaks at higher stimulus frequencies (50-100 ms inter-stimulus interval; Figure 2H) without the need of signal deconvolution. Moreover, neither the amplitude nor the kinetics of the iGluSnFR responses were affected by the competitive AMPA receptor antagonist NBQX (n = 4, Figure 2I), which may offer some flexibility to avoid specific experimental complications, such as minimizing excitability for experiments in highly recurrent networks or minimizing plasticity induction by repetitive and prolonged stimulation paradigms. Altogether, the iGluSnFR method for quantal analysis offers more experimental flexibility and faster kinetics than that afforded by NMDAR-mediated calcium influx detected by GCaMP6f.

2.4. Biophysics of Glutamate Release Variability

The goal of quantal analysis is to infer release properties of glutamate release from a distribution of recorded release magnitudes. Quantal analysis of synaptic release has been performed for decades and the formalism has evolved and adapted as new and improved recording technologies were developed. For didactic purposes, we revisited here some of the basic assumptions commonly held for performing quantal analysis of glutamate release events at single synapses.

We started by exploring the most appropriate continuous distribution to describe the inherent variability expected of a glutamate quantum. Our aim was to derive from the biophysical features of synaptic vesicles a mathematical description of the expected distribution of glutamate release amplitudes, along with their expected variability. Based on previous theoretical studies, we expect that the variability of inner vesicular volumes (Figure 3A) will be a potent determinant of the variability in the amount of glutamate molecules per quantum (Bekkers et al., 1990). What variability of glutamate release do we expect from fluctuations in vesicle diameters only? In order to find this, we first constrain the concentration of glutamate within synaptic vesicles to be constant across the many synaptic vesicles of a given neuron. In addition, we assume faithful release of a single vesicle (in this case release probability p = 1) and that the relative location and loading of vesicles does not introduce a significant amount of variability in the activation of post-synaptic receptors. We will revisit each of these assumptions sequentially as we assemble the mathematical synapse model.

Electron microscopy studies have shown that the variability in vesicle diameters at hippocampal synapses is normally distributed. Using the measurements obtained from one such study (mean vesicle diameter 38.7 nm, CV_d = 0.13 Qu et al., 2009) we generated a simulated distribution of 10,000 inner vesicle diameters (Figure 3B) and a corresponding distribution of the inner vesicle volumes (Figure 3C), assuming the shape of synaptic vesicles is approximated by a sphere. Inner vesicle diameters and volumes were calculated by first subtracting the thickness of the plasma membrane (12 nm Qu et al., 2009). This volume distribution can be readily calculated by a change in variable of the diameter distribution (Bekkers et al., 1990; Barri et al., 2016). In line with the cubic relationship between volume and diameter, the resulting distribution (Figure 3C) is non-Gaussian as it displays an important rightward skew.

To compare the possible distributions of volumes emanating from a range of experimentally derived vesicular diameter, we explored a set of continuous distributions (normal, gamma, Weibull, lognormal) that could accurately describe the skewed distribution of inner vesicle volumes simulated. Using the Bayesian Information Criteria (BIC) as a scoring metric, we ranked the distributions with the degree with which they fit the simulated distribution (Figure 3D). We found that the gamma distribution provided the best approximation of the empirical distribution of vesicle volumes, followed by the Weibull, lognormal, and finally the normal distribution. By assuming equal loading of vesicular glutamate concentration into vesicle, we would therefore expect that the distribution of glutamate content per vesicle would be dictated by that of vesicle volumes (that is, a gamma distribution). These findings are intriguing when we consider that many previous studies of quantal analysis have reported using a Gaussian mixture model of release events (Larkman et al., 1997; Hardingham et al., 2010; Malagon et al., 2016; Jensen et al., 2019), although some have used skewed distributions (Lavoie et al., 2011; Barri et al., 2016) and at least one study a gamma distribution (Bhumbra and Beato, 2013).

The gamma distribution is described by two parameters: a shape parameter γ and a scale parameter λ and it is expressed in terms of the gamma function Γ(·). When used to approximate the distribution of vesicle volumes, v arising from normally distributed diameters, we write

where λ_v and γ_v are the parameters for the volume distribution. These paremeter values can be found by matching the first two moments of simulated (Figure 3) and theoretical (Equation 1) distributions. Equation 1 has a right-skew controlled by the parameter γ_v. Conveniently, its mean (E[v] = γ_vλ_v), its variance ($\mathrm{\text{Var}}\left[v\right]={\gamma }_v{\lambda }_v^2$), its skewness (skewness = $2/\sqrt{{\gamma }_v}$) and its coefficient of variation (CV_v = $1/\sqrt{{\gamma }_v}$) are simple expressions of these parameters. Also of considerable practicality, the addition of two independent gamma-distributed variables results in a random variable that is itself gamma-distributed with shape parameter equal to the sum of the shape parameters. As pointed out by Bhumbra and Beato (2013), these properties allow for a clearer treatment glutatmate release variability without explicitly compromising the validity of the gamma-release model.

To relate the parameters of the gamma distribution with vesicle dimensions, we calculated the expected range of γ_v and λ_v as a function of the mean vesicle diameter (μ_d, Figure 3E) and diameter coefficient of variation (CV_d; Figure 3F) for simulated vesicle volume distributions. The shape parameter is unaffected by changes in mean diameter, but the scale parameter increases nonlinearly with increasing diameters. In addition, the shape parameter decreases and the scale parameter increases when the CV of vesicle diameter increases. It is therefore possible to interpret an increase of the scale parameter as an increase in the mean vesicle diameter, but only if the shape parameter shows no concomitant changes.

What are the theoretical predictions of vesicle volume variability for optical measurements of cleft glutamate? Using the mean diameter μ_d and the variability of diameters CV_d from electron microscopy recordings, we predict λ_v = 0.15 and γ_v = 6.8. Importantly, these parameters give rise to a variability of volumes CV_v of 0.38. In theory, unequal loading of vesicular glutamate content, neurotransmitter diffusion and observational noise should increase the coefficient of variation once we consider the glutamate reported on the post-synaptic membrane instead of vesicle volumes. Since these factors are likely to be captured by another skewed distribution (Franks et al., 2003; Bhumbra and Beato, 2013; Bird et al., 2016), it is appropriate to use a gamma distribution to capture the total variability of univesicular releases. To consider a possible discrepancy between the variability of univesicular releases and that of vesicle volumes, we use γ and λ to parameterize the distribution of univesicular releases, not to be confused with the parameters of the theoretical volume distribution γ_v and λ_v. In fact, since additional sources of variability can only increase the CV, our γ_v should be considered an upper bound on γ. To summarize, we derived biophysical constraints for the parameters of a gamma-distributed set of univesicular glutamate release events (UVR) using an experimentally-derived distribution of inner vesicular volumes and the assumption of equal glutamate loading across vesicles of different sizes.

2.5. Observational Error and iGluSnFR Transduction

In principle, the experimental readout expected from the non-uniform distribution of cleft glutamate will arise in part from the cubic transform outlined above but it can be corrupted by loading, diffusion and by non-optimality of the iGluSnFR signal transform. In order to begin addressing the issue of iGluSnFR transform, we sought to experimentally interrogate as directly as possible the relationship between the quantity of glutamate release at single spines and the amplitude of iGluSnFR-mediated transients. By varying the amount of glutamate released onto dendritic spines through step-wise increments in uncaging laser power during simultaneous optical and electrophysiological recordings (Figure 4A), we found that the relationship between uncaging laser power and iGluSnFR amplitude was linear (Figure 4B) within the expected physiological range of glutamate release, as determined by the average amplitude of uEPSCs (Bé¨ıque et al., 2006; Soares et al., 2014, 2017; Lee et al., 2016). Furthermore, we found an inverse relationship between the uncaging stimulus intensity and the CV of IGluSnFR responses (Figure 4B). Taken together, iGluSnFR can linearly report changes in glutamate concentration at dendritic spines with high precision.

We then estimated a convolved metric of observational error CV_opt to be 0.15, by measuring the variability of the iGluSnFR transients upon presentation of nominally fixed amounts of glutamate concentrations by repetitive uncaging at a fixed laser intensity (around 30 mW; Figure 2E); while uEPSC amplitudes were within an expected physiological range (Figure 4E). At most, adding this measurement noise brings the combination of diameter and optical variability to the upper bound $\sqrt{C{V}_v^2+C{V}_{\text{opt}}^2}=$ 0.40. The formalism outlined above therefore predicts the distribution of optical signals when glutamate is released from a presynaptic terminal. We next considered the variability imposed by the stochastic nature of vesicle releases.

2.6. Release Failures

Large amount of variability is attributed to the stochastic failure of vesicle release (Calvin and Stevens, 1968) upon action potential arrival. To formally include this process in our predicted distribution of optical signals, we considered a mixture model wherein we stochastically sampled from one of many independent sub-distributions, which are called components. Since in certain conditions, multiple vesicle release (MVR) occur at central synapses (Wadiche and Jahr, 2001; Oertner et al., 2002; Conti and Lisman, 2003; Christie and Jahr, 2006), we consider a MVR model for which univesicular release (UVR) is a special case. When n vesicles are docked and ready to be released and when each of these vesicles are released independently with probability p, the number of vesicles released will follow the binomial distribution. It is important here to clarify that the parameter p is not to be confused with the synapse's overall release probability (P_r) – the probability that any of the n docked vesicles will release. It is also important to note that the formalism does not specify whether or how multivesicular release is distributed in nanodomains. Irrespective of their sub-micron localization, we will expect that at times all vesicles have failed to release, in which case we will sample from the failure distribution. Assuming a Gaussian measurement noise – here called optical – for the failure distribution, we obtain the gamma-Gaussian mixture for observations of fluorescence amplitude F

where ${\sigma }_{opt}^2$ is the variance of the optical noise (derived in Section 2.5) and G(F|μ, σ²) is a Gaussian distribution of mean μ and variance σ². In Equation 2, k ranges from 1 to n and refers to the possible number of vesicles released. The binomial coefficient p^k(1 − p)^n−k establishes the probability of observing k vesicles, while each time that k vesicles are released, we sample from a gamma distribution g(F|kγ, λ) with a shape parameter corresponding to k times the univesicular shape parameter γ.

We make three observation on this gamma-Gaussian mixture model of glutamate release at single synaptic contacts. Firstly, we distinguish the vesicular release probability p from the probability of any vesicle being released $P_r=1-{(1-p)}^n$. Secondly, the mean and the variance of this distribution now depend on the maximum number of vesicles released n, namely μ = npγλ and ${\sigma }^2={\sigma }_{opt}^2{(1-p)}^n+{\lambda }^2\gamma np(1+\gamma (1-p))$. Lastly, it can be useful to analyze the measured variability, CV, in terms of the variability of univesicular releases $C{V}_{\text{UVR}}=1/\sqrt{\gamma }$, the variability due to observational error σ_opt and the variability of a binomial process $C{V}_{\text{bin}}^2=(1-p)/np$. In this way we can parse the variability in three terms (see Methods for derivation)

This expression allows us to parse out the variability in terms of distinct sources.

Overall, for the experiments described in Figure 1, the gamma-Gaussian mixture should capture the variability of glutamate-dependent optical events originating from: Optical (various optical measurement noise), Binomial (the stochastic behavior of releasing n vesicles independently), and Unitary (release variability of associated with each vesicle release). The latter comprises variability from vesicle sizes, loading and diffusion. It has a total of five parameters: ${\sigma }_{opt}^2$ the variance of the optical noise, n the number of vesicles, p the probability of each vesicle being released, λ the scale and γ the shape parameters of the gamma distribution.

2.7. Inferring Release Parameters From Quantal Peaks

An intuitive approach to discriminate release events from failures lies in classifying a trial as a failure of release if the observed peak fluorescence is less than twice the standard deviation of the optical noise, and success otherwise (Figure 5A). From the distribution of success amplitudes, one then extracts the mean and coefficient of variation, called potency and CV_suc, respectively. It is not immediately clear, however, how false positives and false negatives arising from a thresholded detection method influences the estimates of potency and CV_suc. In this section, we use computer simulations to determine the bias introduced by optical noise on these measures.

To quantify the bias arising from classification errors, we generated surrogate amplitudes and calculated the potency and CV_suc using a threshold corresponding to two standard deviations of the optical noise (Figures 5A,B). We compared these estimates to potency and CV_suc calculated without classification errors. We found that for γ = 2 and λ = 0.15, classification errors leads to an over-estimation of the potency for all release probabilities (Figure 5C). This overestimate was restricted to the lower range of shape-parameter (Figure 5D) and scale parameter values (Figure 5E). These biases are overall relatively small, but the effects of optical noise are more dramatic on the calculation of CV_suc. Given γ = 2 and λ = 0.15, we found that CV_suc is drastically underestimated for all p (Figure 5F). This underestimate arises in a range of shape and scale parameters with low values (Figures 5G,H). In the case of threshold classificiation of successes and failures, we conclude that CV_suc will be heavily underestimated when the skewness is noticeable and the quantal size (γλ) is small.

Next we investigated the the consequence of skewed distribution on the identification of quantal parameters n and p. Common approaches to estimate quantal parameters are based on the identification of quantal peaks (Larkman et al., 1991, 1997; Hardingham et al., 2010; Malagon et al., 2016). These approaches assume that the observation of a peak in the release-amplitude histogram can be read off as a quantal mode, an assumption that is often difficult to justify (Clements, 1991; Walmsley, 1995; Ninio, 2007). Peak identification can be even more problematic when the release components show an important skew. Indeed, we noted that mixtures of skewed distributions rarely show quantal peaks (Figure 6). For instance, a gamma-Gaussian mixture with n = 2 will transition from the absence of quantal peaks (Figure 6A) to the presence of quantal peaks (Figure 6B) only if the skewness of the components is reduced beyond the range predicted from biophysics (Figure 6C). These observations extend the limitations previously raised (Clements, 1991; Walmsley, 1995; Ninio, 2007) and show that analysis of quantal peaks is problematic especially when the distribution of univesicular release is skewed or only for a very narrow range of release probability. Since we expect a significant skew from known vesicle diameters (Figure 3), we sought a different method for extracting release properties.

2.8. Inferring Release Properties Using Likelihood Maximization

Maximizing the likelihood function provides an appealing alternative to feature-based approaches such as Bayesian quantal analysis (Bhumbra and Beato, 2013), non-stationary fluctuation analysis (Silver, 2003; Evstratova and Tóth, 2014), or to quantal-peaks approaches (Larkman et al., 1991). This approach does not rely on a trial per trial classification of successes and failures. Instead, the task is to find the set of parameters (n, p, γ, λ, σ_opt) that maximizes the probability of observing all recorded amplitudes given our gamma-Gaussian model. In the case of the likelihood written in Equation 2, there is no guarantee that only one such maximum exists, which means that it can be difficult to find the global maximum. Likelihood maximization algorithms can greatly help in solving this type of task.

For our problem, the Expectation-Maximization (EM) algorithm appears a natural choice since it was developed to improve parameter inference in mixture models (Dempster et al., 1977). The EM algorithm has been used previously to infer synaptic properties, but using different experimental and computational methodologies (Barri et al., 2016). For efficient use of this algorithm, it is critical to derive estimation formulas specific to a given problem. Since we are not aware of any such estimation formulas for the gamma-Gaussian mixture (Equation 2), we next describe our adaptation of the EM algorithm.

The likelihood maximization in the EM algorithm is associated with the principle of gradient ascent (Xu and Jordan, 1996). Accordingly, it begins with an initial guess, and then iteratively updates these estimates to gradually maximize the likelihood $L(\text{F}|\theta ,n,{\sigma }_{\text{opt}}^2)$ of observing the N observations of fluorescence amplitude denoted by the vector F given the parameters θ = (p, γ, λ). Using an initial guesses θ₀ = (p₀, γ₀, λ₀), the algorithm will have found the optimal value of each parameter $\widehat{\theta }=(\widehat{p},\widehat{\gamma },\widehat{\lambda })$. The parameter n will be treated as a meta-parameter to the EM algorithm, whose optimum is obtained by finding the $\widehat{n}$ with its own optimal $\widehat{\theta }$ that maximizes the likelihood $L(\text{F}|\widehat{\theta },n,{\sigma }_0^2)$, or equivalently, minimizes the negative log-likelihood. The variance of the optical noise, ${\sigma }_0^2$, can be estimated independently by calculating the variance of the null distribution (see Methods).

Typically, a good initial estimate of the parameters can greatly speed the inference process. In the present case, we have argued that a good prior on the shape parameter can be obtained from the biophysics of vesicle release with known, Gaussian distributed, vesicle diameters (Qu et al., 2009). We initialize the shape parameter to a value of γ₀ = 4. To initialize the probability of release, we observe that only optical noise can capture fluorescence amplitudes smaller than zero. Therefore, we compute the fraction, c, of the total number of observations falling below zero and equate this to half the failure probability. This suggests the initialization p₀ = 1 − (2c)^1/n. There remains the initial value of the scale parameter. Given that the mean of fluorescence amplitude of the model is npγλ, we use the mean of the observed fluorescence amplitudes μ_F to initialize λ₀ = μ_F/np₀γ₀.

The EM algorithm is iterative and variational. That is, it first approximates the likelihood by an auxiliary function, which we will call Q. It then iterates between a maximization of this auxiliary function (the maximization step) and an improvement to the approximation by generating a new auxiliary function (the expectation step). Using b(k|n, p) to denote the binomial distribution, the likelihood over N observations

is replaced by

This auxiliary function relies on N(n + 1) variables ${\mu }_{i,k}^{(t)}$. These are the posterior probabilities of sampling from the kth component given a guess of the parameters θ^(t).

In the expectation step, we compute the posterior probabilities for k > 0

The posterior probabilities for k = 0, ${\mu }_{i,0}^{(t)}$, would need to be considered only if we were to use the EM algorithm to determine σ_opt. Importantly, these posterior probabilities are computed using the previous guess θ^(t) = (γ^(t), λ^(t), p^(t)).

In the maximization step, we compute the new parameter set, which maximizes the auxiliary function ${\theta }^{(t+1)}={\mathrm{\text{argmax}}}_{\theta }Q(\theta ,{\theta }^{(t)})$. This is done via three re-evaluation formulas, obtained by setting the gradient of Q to zero. In what follows, we will use ${H}(·)$ to denote the Heaviside function. The first formula gives an update of p

To compute the second, we first calculate the model mean

and then maximize the terms of Q that depend on γ^(t+1)

The third formula updates the scale factor

The expectation and maximization step are then repeated in alternation until convergence, which is defined by a tolerance value on the likelihood update.

We use these parameter estimates to compute the log-likelihood using Equation 4. Repeating the EM-method for n within a physiological range of 1-10 allows us to find the number of vesicles $\widehat{n}$ which maximizes the log-likelihood

Since the results may depend on the initialization point, we repeat the procedure with ten different initialization points. The parameter values associated with the highest likelihood become our parameter estimates.

2.9. Statistical Inference on Surrogate Data

To determine the precision and the validity of the EM method for extracting release properties, we apply the method on simulated data. We assume that the fluorescence amplitude are sampled from the gamma-Gaussian distribution. Once a sample is drawn, we will use the EM method to extract the release properties, namely the parameters γ, λ, n and p. Knowing the true parameters allows us to calculate the average difference between estimated and true parameters (bias) and the size of random fluctuations in the estimated parameters (variance). Since these estimator bias and variance will depend on the specific set of parameter values used to generate surrogate data, we must explore different types of parameter values. For the sake of illustration, we consider three cases: i) Univesicular release (Figure 7A), ii) multivesicular release with a low value of the shape parameter corresponding to the absence of dip in the probability distribution (Figure 7B) and; iii) multivesicular release with a high value of the shape parameter leading to well resolved quantal peaks but inconsistent with the biophysical constraints (Figure 7C).

We computed the bias and standard deviation of the estimates using 500 surrogate experiments and the expectation-maximization algorithm of the gamma-Gaussian mixture (see Methods Sect. 4.6). Since both the bias and the standard deviation are expected to depend on the number of samples per dataset, N, we report the bias and standard deviation as a function sample size. The correlation coefficients shown in Figure 7D reveal two interactions. Firstly, release probability as well as the shape and scale parameter estimates are strongly correlated. Secondly, these three parameter estimates are anti-correlated with the estimate of the number of vesicles. These compensations are also reflected in the sample-size dependent biases, where an underestimate (overestimate) in n is accompanied by an overestimate (underestimate) in the other parameters (Figures 7E–H). This reflects the fact that n is determined in a separate step from the other parameters and that for a larger n the other parameters must decrease to keep the same mean amplitude. We have verified with simulated sample sizes of 50 000 that these biases are restricted to small sample sizes, at very large sample sizes the biases reach zero.

Next, we consider the bias and variance of estimators for a sample size of 50, which represents a realistic sample size for our experimental conditions. At those sample sizes, we find that our method underestimates release probability, shape parameter and scale parameter (Figures 7E–G, red lines). These biases reflect the fact that, given the small number of vesicles considered, the number of vesicles can only be overestimated (Figure 7H). Considering sample size of 50 but for surrogate data with two vesicles, there remains a small underestimation of the shape parameter (Figures 7E–G, green lines) but the bias in the number of vesicles is much reduced (<0.25), and so is the bias in release probability and scale parameter. These biases are further reduced and become negligible in the less realistic situation where quantal peaks can be identified (Figures 7E–H, blue lines). Lastly, we note that the estimator standard deviations at N = 50 are sufficiently substantial to require averaging over multiple synapses in order to make precise parameter estimates.

2.10. Statistical Inference on Experimental Data

We next apply this EM algorithm on experimental data from iGluSnFR-mediated optical recordings of glutamate release. We used recordings of iGluSnFR-mediated signal induced by trains of ten axonal electrical stimulation at low frequencies (1,2,4, and 8Hz), from which we extracted a distribution of fitted release amplitudes (see Methods). Amplitude distribution from an exemplar spine is shown in Figure 8A. This distribution is captured very well by the gamma-Gaussian mixture model (Equation 2). The best fit for this recording was achieved for shape parameter γ = 1.4, scale parameter λ = 0.2, release probability p = 0.42 and 2 vesicles. Figure 8A shows that the theoretical distribution fits the empirical histogram well. This fit arises from individual components having an important skew. The inset of Figure 8A shows the negative log-likelihood as a function of the number of vesicles. Although there is a clear minimum at n = 2 vesicles, the curvature is fairly large, as is predicted by the small estimator variance (Figure 7H) under Cramer-Rao inequality. Importantly, the likelihood is considerably worse for the n = 1 model compared to any n > 1 models. Altogether, parameter inference using this EM algorithm on iGluSnFR-based analysis of glutamate release on single synapse argues that an action potential stochastically triggers the fusion of a few vesicles releasing a variable and highly skewed amount of glutamate.

We repeated this analysis on a set of experiments from 18 different spines (Figure 8B). Here, the average number of vesicles fitted by the algorithm was 2.10 ± 0.3, while individual spines were best fit by n ranging between 1 and 3 vesicles. The average shape parameter value was 1.42 ± 0.08. This parameter regulates the univesicular releases, the univesicular CV_UVR, described in the biophysical predictions. In principle, the average shape parameter fitted by the EM algorithm should correspond to a univesicular CV_UVR of 0.84, but we recall that our estimates of the shape parameter were shown to bear a small-sample bias, which we estimated to negative 0.2 (Figure 7F). As a consequence, our bias-corrected estimate of univesicular CV is 0.77. In comparison, we predicted that a CV_v = 0.38 (Figure 3) would arise from known vesicle dimensions, thus a difference of 0.39. To see how our estimate of univesicular CV depends on the number of vesicles in the model, we fixed n and inferred CV_UVR for each spine (Figure 8C). We find that increasing n increases the CV_UVR inferred. This CV_UVR remained high and above both the variability expected from volumes (CV_v) and the variability of thresholded successes (CV_suc). This is consistent with the view that CV_v is a lower bound (Lavoie et al., 2011) and CV_suc is underestimated (Figure 5). In sum, statistical inference of our gamma-Gaussian model suggests CA3-CA1 stochastically release 1-3 vesicles with variable quantum.

The formalism outlined above allows to begin parsing out the variability of synaptic transmission at single synapses. Using the average parameters extracted using the expectation-maximization algorithm, Equation 3 can be used to separate the variability of observed evoked amplitudes in three terms. The first term captures the variability due to the glutamate sensor itself and to concurrent optical measurements The second term captures the fluctuations in release amplitude attributable to a single vesicle release and scaled by the average number of vesicle released (i.e., diffusion, loading and vesicle volume). The last term captures the variability of releasing sometimes two, one or zero vesicles (with zero unitary variability). We named these sources of variability optical, unitary and binomial, respectively. As shown in Figure 8D, we estimate that that 4% of the variability was optical, 30% binomial and 66% was unitary. Thus the results suggests that, despite the fluctuating number of vesicles being released, the variability of synaptic transmission arises mostly from the variability in unitary vesicle content released.

3. Discussion

The use of the glutamate fluorescent reporter iGluSnFR provides a valuable proxy of glutamate release at single central synapses. To interpret the variability of glutamate release observed in our recordings, we built a gamma-Gaussian mixture model based on stochastic release of vesicles, each with a variable diameter size and additional sources of variability. We highlighted important biases in the traditional measures of the variability of successes and provided an alternative method based on the expectation-maximization algorithm. Our statistical method showed small biases on surrogate data and allows inference of estimates of quantal parameters. Together, these experimental and analytical tools allow to resolve the dynamic structure of synaptic transmission.

Optical quantal analysis confers several distinct advantages over classical electrophysiological methods, but also some limitations. One of the main advantages of optical methods is that the experiment is localized at an unambiguous source spine, thereby removing classical uncertainties such as the unknown location of synaptic inputs, the impact of dendritic filtering, and the determination that a single synaptic input is being monitored. While strong criteria have been developed in the past to classify electrophysiological data as arising from a single synaptic contact (i.e., minimal stimulation criteria; Malinow and Tsien, 1990; Raastad et al., 1992; Stevens and Wang, 1995; Dobrunz and Stevens, 1997, 1999), some of these criteria may actually introduce false rejections and selection biases in population sampling, favoring against synapses that display multi-quantal, highly-variable, and/or a high probability of release. Optical methods are not without their disadvantages, chief amongst them being their still poorer temporal resolution as compared to traditional electrophysiological methods. Nevertheless, advances in molecular engineering have enabled faster and more specific optical reporters of synaptic transmission, including iGluSnFR which offers a significant temporal improvement over some of the fastest genetically-encoded optical calcium reporters (Figure 2E). Interestingly, the iGluSnFR family of reporters is still growing to include faster variants and variants with different emission spectra (Marvin et al., 2018; Wu et al., 2018).

One feature of the iGluSnFR relies on its non-reliance on postsynaptic glutamate receptor activation. Interpretation of calcium-based and electrophysiological-based measurements of synaptic release that rely on glutamate receptor activation are confounded by issues such as the non-linear relationship between glutamate concentration, glutamate receptor conductance and calcium-mediated fluorescence (Smith and Howe, 2000), the non-uniform distribution of glutamate receptors in the postsynaptic membrane (Biederer et al., 2017), and the distance between the presynaptic site of release and postsynaptic receptors (Franks et al., 2003), all of which are difficult to measure and contribute to convolving the end-result signal. iGluSnFR, being plasma membrane localized but lacking postsynaptic anchoring domains, is presumably evenly distributed on the cell surface and reports glutamate release largely independently of the precise location of vesicle release. Additionally, iGluSnFR reports glutamate release events in the presence of glutamate receptor antagonists (Figure 2I) which offers experimental flexibility. Finally, as demonstrated in Figure 3F, iGluSnFR provides a linear report of physiologically-relevant glutamate release, which allows for a more direct quantitative interpretation of the glutamate signal. Altogether, the experimental method we describe addresses a number of historical limitations and provides a welcome complement to existing methodologies to study basic features of synaptic transmission and plasticity.

Synaptic transmission is variable. Obtaining an accurate estimate of the size of this variability is an obligatory step in order to parse information content from noise during neural communication. By applying a traditional threshold-based classification of successes and failures on iGluSnFR transients, we obtained a fairly low average variability CV_suc ≈0.5 (Figure 8C). Some of our recordings showed CV_suc in the 0.2-0.4 range, which closely matches the values reported for putative single-synapse electrophysiological recordings using either manual or threshold-based classification of failures and successes (Bekkers and Stevens, 1995; Dobrunz and Stevens, 1997; Hanse and Gustafsson, 2001). We however readily observed synapses that showed a higher CV_suc (up to 0.8) when optically probed. It is likely that these synapses would have been ignored when applying selection criteria commonly used for minimal stimulation experiments. Moreover, threshold-based classification inherently introduces classification errors, which can dramatically alter estimates of CV_suc. The statistical methodology presented here should circumvent this issue. Consistent with our estimates on surrogate data (Figure 5), we observed that the variability of individual synaptic release can be much higher than threshold-based CV_suc and reach CV_UVR = 0.8. Further consistent with the effect of classification errors, our estimate of the average release probability of individual release is higher using the expectation algorithm (p = 0.69; Figure 8B) than using a threshold-based approaches (previous estimates were <0.61 Hanse and Gustafsson, 2001, 0.4 Bekkers and Stevens, 1995 and 0.2–0.4 Dobrunz and Stevens, 1997). Our assumption that glutamate is packaged at equal concentration in vesicles may be factually incorrect, implying that neurons have evolved a means to normalize glutamate content across vesicles that differ in volume. Such a result imply a redistribution of the relative breakdown of factors that contribute to quantal variability Figure 8D. Irrespective of this interpretation, our results further suggest that vesicle fusion can release very small but non-zero glutamate transients and that glutamate release variability is considerably higher than previously thought.

The usefulness of the experimental methodology and analytical formalism described herein extends beyond our preliminary efforts to identify the source of quantal variability. The iGluSnFR family of optical reporters are well suited to study the dynamical mechanisms that regulates synapse-specificity such as those controlling glutamate spillover (Asztely et al., 1997; Chalifoux and Carter, 2011; Lee, 2012) and its downstream impacts on neuronal function. Furthermore, by providing a direct proxy of release, it is well suited to interrogate several features of synaptic plasticity mechanisms across a wide range of central synapses and experimental preparations, including in vivo. In principle, optical fluorescent reporters may also be amenable to use in combination with complementary imaging modalities: for instance, one can envision using super-resolution imaging to estimate how strongly the modeled parameter n – the maximum number of vesicles released simultaneously – relates to optical estimates of the total number of active release sites (i.e., transsynaptic nanocolumns). Such investigations could provide a deeper understanding of the precise nature of quantal variability and its impact on the plasticity of information transfer at central synapses.

4. Methods

The essential elements of optical quantal analysis are described in the main text. In this section, we give additional precisions on experimental, analytical and computational methods.

4.1. Organotypic Slices and Biolistic Transfection

A detailed description of our methodology for hippocamal organotypic slice preparation and biolistic transfection is described in (Soares et al., 2014). Briefly, organotypic slices were prepared from both male and female postnatal day 7 Sprague Daley rats (Charles River Laboratories, MA, USA) using the interface method originally described in Stoppini et al. (1991). In accordance with protocols approved by the University of Ottawa's Animal Care Committee, animals were anesthetized in an isofluorane infused chamber, decapitated, and hippocampi were removed in ice cold cutting solution containing (in mM): 119 choline chloride, 2.5 KCl, 4.3 MgSO₄, 1.0 CaCl₂, 1.0 NaH₂PO₄–H₂O, 1.3 Na-ascorbate, 11 glucose, 1 kynurenic acid, 26.2 NaHCO₃, saturated with 95% O₂ and 5% CO₂ (pH = 7.3; 295-310 mOsm/L). Transverse slices were cut at 400 μm thickness using a MX-TS tissue slicer (Siskiyou, Grants Pass, OR) and cultured on 0.4 mm millicell culture inserts (EMD Millipore, Etobicoke, Canada) at a temperature of 37^oC.

Hippocampal slices were transfected at 6-7 DIV using a hand held gene gun (Biorad, Hercules, CA). Cartridges for the gene gun were prepared in advance by precipitating 50 μg of cDNA plasmid onto 8-10 mg of gold microparticles (1.0 μm diameter; Biorad) at a ratio of 80/20 by weight of either iGluSnFR or GCaMP6f and mCherry cDNA plasmid, respectively. The precipitation step was performed in a 0.1 M KH₂PO₄ buffer solution containing 0.05 mM protamine sulfate (rather than spermine, as per previous protocols). The DNA-gold precipitate was washed and suspended (3 times) in 100% ethanol before loading in the tubing station. Once the cartridges were dried and cut, they were placed in a sealed container with desiccant pellets at 4 ^oC until used. The DNA-coated gold particles were delivered to the slice using ~180 psi of helium air pressure. A modified gene gun barrel was used to protect slices from helium blast (Soares et al., 2014). Imaging experiments were performed 3-5 days after biolistic transfection.

4.2. Optical Recording of IGluSnFR Transients

Slices were removed from culture and placed in a custom recording chamber under a BX61WU upright microscope (60X, 1.0 NA objective; Olympus, Melville, NY). Slices were continuously perfused with a Ringer's solution containing (in mM): 119 NaCl, 2.5 KCl, 4 MgSO₄-7H₂0, 4 CaCl₂, 1.0 NaH₂PO₄, 11 glucose, 26.2 NaHCO₃ and 1 Na-Ascorbate, saturated with 95% O₂ and 5% CO₂ (295-310 mOsm/L). For evoked stimulation experiments, a glass monopolar electrode filled with Ringer's solution was positioned adjacent to transfected cells in the direction of CA3. Simultaneous two-photon imaging of iGluSnFR and mCherry was performed using a Ti:Sapphire pulsed laser (MaiTai-DeepSee; Spectra Physics, Santa Clara, CA) tuned to 950 nm. Emission photons were spectrally separated using a dichroic mirror (570 nm) and the emitted light was additionally filtered using two separate bandpass filters (iGluSnFR: 495-540; mCherry: 575-630). The sampling frequency of our line-scan experiments depended on length of the imaged line segment (drawn over a spine and its parent dendrite), but was typically in the range of 1.2 - 1.5 ms / line for all optical quantal analysis experiments. This sampling rate was more than sufficient to fully capture and quantify the rise and falling phases of iGluSnFR transients. In our hands, an optimal trade-off between signal-to-noise, sampling frequency, and reduced bleaching, was obtained by using a 4 μs pixel dwell time. In the frame scan configuration, the sampling limit of our optical system was 65 ms/frame (2 μs pixel dwell time; 256 x 256 pixel frame) when scanning bidirectionally, which was sub-optimal for optical quantal analysis.

Surveying methodology was designed to increase the probability of finding dendritic spines that were responsive to the electrical stimulus. While our frame scan configuration offered the spatial resolution to monitor several spines at once, we found it difficult in practice to identify rapid iGluSnFR-mediated transients due to a low signal to noise ratio and low sampling rate. As a result, line scans were exclusively used to survey the dendritic arbor for responsive spines. Short duration (0.1 ms) low intensity (5-25 mA) stimuli were delivered to the slice at low-frequency (0.1 Hz) while randomly surveying dendritic spines in the apical arbor of transfected cells. To facilitate the process of finding a responsive spine, line scans were performed simultaneously through multiple nearby dendritic spines and, generally, a paired-pulse stimulus (50-100 ms inter stimulus interval) was delivered to increase the probability of detecting glutamate release. Dendritic spines that were unresponsive to an initial probing phase consisting of 5-10 paired pulse stimuli, were not considered for further analysis, while spines demonstrating responsiveness to these initial probing stimuli were selected for quantal analysis experiments. Fluorescent transients were resolvable by eye and on-line analysis was not necessary. Prior to starting an optical quantal analysis experiment at a responsive spine, the stimulus intensity was gradually reduced up to a minimum where time-locked responsiveness was still observed. The process of identifying a responsive spine was generally not trivial and often necessitated several re-positioning of the stimulating electrode. Once a responsive spine was found, however, it was extremely rare to lose fluorescent responsiveness in response to electrical stimulation during an experiment.

4.3. Whole-Cell Electrophysiology and Two-Photon Glutamate Uncaging

Whole-cell recordings were carried out using an Axon Multiclamp 700B amplifier. Electrical signals were sampled at 10 kHz, filtered at 2 kHz, and digitized using an Axon Digidata 1440A digitizer (Molecular Devices, Sunnyvale, CA). Transfected CA1 pyramidal neurons were targeted and patched using borosilicate glass recording electrodes (World Precision Instruments, Sarasota, FL) with resistances ranging from 3-5 MΩ. All uncaging evoked currents were recorded at a holding potential of −70 mV. The intracellular recording solution contained (in mM): 115 cesium methane-sulfonate, 0.4 EGTA, 5 tetraethylammonium-chloride, 6.67 NaCl, 20 HEPES, 4 ATP-Mg, 0.5 GTP, 10 Na-phosphocreatine (all purchased from Life Technologies, Carlsbad, CA) and 5 QX-314 purchased from Abcam (pH = 7.2-7.3; 280-290 mOsm/L). The extracellular solution was similar to the one described above but also contained 2 mM MNI-glutamate-trifluoroacetate (Femtonics, Budapest, Hungary) and a reduced concentration of MgSO₄-7H₂0 (1.3 mM) and CaCl₂ (2.5 mM). For glutamate uncaging experiments, a second laser line tuned to 720 nm was used to deliver 1 ms light pulses to the tips of dendritic spines while the other laser was tuned to 950 nm to image the uncaging-evoked iGluSnFR transients.

4.4. Analytical Derivation for Three Terms of Variability

We consider the gamma-Gaussian model described by Equation 2 of main article. To determine CV² we note that the mean of Equation 2 is μ_x = npγλ, the mean of the number of vesicles is $\sum _{k=0}^{n}kb(k|n,p)=np$ and the mean amplitude per vesicle is ∫xg(x|γ, λ)dx = γλ. The variance of x is

The first term of this equation becomes

and the second term is evaluated by centering the quadratic part of the integrand around the mean of a component, kγλ. This gives for the second term of Equation 12

To evaluate this expression, one can isolate the variance of the number of vesicles: $\sum _{k=0}^n{(k-np)}^{2}b(k|n,p)=np(1-p)$, as well as the variance of amplitude from k vesicles ∫ (x − kγλ)²g(x|kγ, λ)dx = γλ². Since the last term vanishes in Equation 14, we obtain

Using ${\text{CV}}_{UVR}^2=1/\gamma$ and ${\text{CV}}_{bin}^2=(1-p)/np$, taking the ratio ${\text{CV}}^2={\sigma }_x^2/{\mu }_x^2$ we obtain Equation 3.

4.5. Regression for Amplitude Extraction

We describe the use of a template to extract the amplitude of the evoked responses. The method involves two steps. First we extract a template time-course k by computing the trial-averaged fluorescence response that is triggered by the electrical stimulation. This template is discretized, starts at the stimulation time and ends at a pre-defined time T after it. Trial averaging is performed on responses sufficiently isolated in time to be exempt from other synaptic events. For each trial i the template is scaled by β chosen so as to minimize the mean-squared error with the observed fluorescence ${\text{F}}_{0:T}^{(i)}$ in the corresponding time window indicated by the subscript 0:T. The solution of this least-square problem is well known and follows

In order to report the maximum of the evoked waveform, we scale the β_i by the maximum value of the template. This is the value reported in Figure 8. We calculated the null distribution by fitting the template on fluorescence measurements without electric stimulation. The variance of the null distribution serves as our estimate of measurement noise ${\sigma }_{opt}^2=0.07$.

4.6. Surrogate Data Analysis

To generate surrogate data, we simulate n surrogate fluorescence amplitudes and estimate the parameter values. Each surrogate experiment is repeated M times in order to have M sets of parameter estimates ${\widehat{\theta }}_j$. Using this set of surrogate experiments, we can compute the bias, the variance and the correlation coefficients of the estimates. The bias is calculated by averaging the difference between the estimated parameter and the simulated parameters

In this way, a bias greater than zero means that the parameter tends to be overestimated, and a bias smaller than zero means that the inferred parameters are erroneously small.

To estimate the precision of the estimates, we calculate the variance:

A small variance means that the estimate is precise, although it may or may not be valid.

In addition, we compute the correlation coefficient between different parameter estimates. The correlation coefficient between a parameter θ^(r) and θ^(q) is simply

where σ_r and σ_q are the square root of the estimator variances for θ^(r) and θ^(q), respectively. When the correlation coefficient is positive, it means that estimation errors tend to covary with the same sign, and when the correlation coefficient is negative, it means that an estimation error in parameter r is associated with an estimation error of opposite sign in parameter q.

5. Additional Requirements

For additional requirements for specific article types and further information please refer to Author Guidelines.

Data Availability

The datasets generated for this study are available on request to the corresponding author.

Author Contributions

CS, J-CB, AL, and RN conceived the study. CS performed the experiments. CS, DT, and RN analyzed the data and performed the simulations.

Funding

This work was supported by CIHR grant 14242 (J-CB) and NSERC Discovery grant 06872 (RN) and 05830 (J-CB) and Brain Canada Support (AL and J-CB).

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.

References