| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
INNOVATIVE METHODOLOGY
Department of Statistics and Center for the Neural Basis of Cognition, Carnegie Mellon University, Pittsburgh, Pennsylvania
Submitted 25 June 2004; accepted in final form 23 March 2005
 |
ABSTRACT |
|---|
|
|
|
INTRODUCTION |
|---|
|
; Johnson 1996; Kass and Ventura 2001, and the references therein). In many contexts, however, the conditions of the experiment, or the internal state of the subject, may vary across repeated trials enough to produce discernible trial-to-trial spike train variation beyond that predicted by Poisson or other point processes. Such trial-to-trial variation may be of interest not only for its physiological significance (Azouz and Gray 1999; Hanes and Schall 1996) but also because of its effects on statistical procedures. In particular, as observed by Brody (1999a,b), Ben-Shaul et al. (2001), and Grün et al. (2003), trial-to-trial variation can affect the assessment of correlated firing in a pair of simultaneously recorded neurons. In this paper we present a statistical procedure for testing and estimating trial-to-trial variation in time-varying firing rate, and apply it to the problem of assessing time-varying dependency between spike trains from two neurons by extending the method of Ventura et al. (2005a).
One aspect of trial-to-trial variation is the tendency for neuronal response to shift in time. That is, a neuron may tend to fire earlier or later on some trials than on others, so that realignment of trials becomes desirable (Baker and Gerstein 2001; Ventura 2004; Woody 1967). It is useful to distinguish such latency effects from variation in the amplitude of firing rate, which is sometimes called "excitability." It is possible to describe trial-to-trial amplitude variation by applying a kernel smoother to each trial's spike train (using a Gaussian filter or something similar; Nawrot et al. 1999) or by smoothing the interspike intervals of each trial and inverting to get a firing rate (Pauluis and Baker 2000). However, smoothing each trial ignores a special feature of this situation: although there is substantial information on the general shape of the firing-rate function obtained from the aggregated trials (i.e., the smoothed peristimulus time histogram [PSTH]), there is relatively little information per trial from which to estimate a time-varying function, at least for smaller firing rates (e.g., <40 Hz). This creates a compelling need for statistically efficient procedures; otherwise, standard errors and confidence intervals become very wide and statistical tests have little power (Kass et al. 2005). To improve estimation and inference, simple representations that take account of the aggregate pattern are needed. One such statistical model was discussed by Brody (1999a,b), who took the firing rate on each trial to be the sum of two components: a background constant firing rate multiplied by a trial-dependent coefficient and a stimulus-induced time-varying function multiplied by a second trial-dependent coefficient. Although likely to capture some dominant features of trial-to-trial variation, this model may be too restrictive for many situations: it requires an experimental period during which the neuron fires at a constant background rate, and it assumes a single multiplicative constant describes the fluctuation in stimulus-induced firing rate. Additional statistical issues concerned the identification of the end of the background period and beginning of the stimulus-induced period, the use of the raw PSTH rather than a smoothed version of it to estimate the stimulus-induced firing rate, and the assumption of Poisson spiking. The method described here avoids the requirement of a background period because it fits the excitability effects as continuous functions of time, allows two or more coefficients to describe the fluctuation in firing rate, and may be applied with non-Poisson spiking. It also may be used to assess whether excitability effects are shared by two or more neurons.
|
|
METHODS |
|---|
|
r(t) of Pr(t) for r = 8, with 95% probability bands, obtained by 1) smoothing (filtering) the data for trial 8 to produce an initial estimate 
r(t), 2) smoothing the PSTH to obtain an estimate (t), 3) smoothing further the residual r(t) – (t) to get the final estimate r(t) as the sum of (t) from step 2 and the smoothed residual from step 3. This procedure is explained in greater detail in the next subsection. We used splines to do the smoothing, allowing a small number of degrees of freedom (large bandwidth) for step 3 but a larger number of degrees of freedom (smaller bandwidth) for steps 1 and 2.
|
four times more efficient than the estimates in D, or equivalently, that the method in D would require, if it were possible, that we use four repeats of the same trial r = 8 to achieve the same accuracy as the estimates in B.
In Fig. 1A, the 20 spiking probability functions Pr(t) are related to P(t) according to
 | (1) |
r is a latency that varies across trials. Thus the variability of the spike trains is attributed partly to the stochastic variation of point processes and partly to the extra variation caused by the gain factor gr and shift (latency) r, which might represent trial-dependent changes in the internal state of the subject. We have introduced the shift parameters r for completeness. We do not discuss here the estimates of latency effects because of the availability of methods for estimating these parameters; they may be obtained, for example, by the procedure of Ventura (2004) as a preliminary step. In the next subsection we address the statistical problem of recovering the factors gr(t). Because the spike train for a single trial is relatively sparse, it is important to reduce the dimensionality of the functional representation. To achieve that, we introduce a principal-component decomposition of the functional variation gr(t), which serves in part to reduce dimensionality and in part to provide additional interpretation of the variability.
We show how significance tests may be used to examine whether a constant-excitability model is appropriate, and whether any additional variability may be summarized by one or two parameters, in terms of the principal components of the variability.
One purpose of analyzing trial-to-trial variation is to examine its effect on time-varying dependency between two neurons. In the remaining subsections we extend the methods of Ventura et al. (2005a) to incorporate trial-to-trial variation; we show how to adjust the cross-correlogram and JPSTH for trial-to-trial variation; and we show how significance tests may be used to produce evidence that two neurons have shared input in the form of shared trial-to-trial variation.
Equation 1 is the basis for the methodology developed here. Firing rate functions that depend on t are valid for Poisson processes only. For non-Poisson processes, a conditional intensity P(t – r|Hrt) must replace P(t – r) in Eq. 1, as in Ventura et al. (2005a), where Hrt is the history of trial r up to time t. An example of a simulated non-Poisson neuron is treated in RESULTS, Simulated data: non-Poisson spike trains.
Low-dimensional representation of time-varying excitability effects
If Eq. 1 is believed to apply, so that the firing rate of a neuron is different across trials, an estimate of P(t), generally the PSTH, will not provide a complete description of the firing properties of the neuron. One standard way to proceed is to fit, by Gaussian filtering or by other means, the firing-rate functions on all trials separately. Because one trial has few spikes, at least relative to the number of spikes available to smooth a PSTH, smoothing each trial separately will produce highly variable estimates. We propose a more efficient alternative.
We first obtain a smoothed estimate (t) of P(t) by smoothing the PSTH. We prefer a spline-based method called BARS (Bayesian adaptive regression splines) (DiMatteo et al. 2001; Kass et al. 2003) for this purpose, but many alternative smoothing techniques would work well for typical data sets. The advantage of BARS is that it determines automatically the optimal number and placement of knots. An analogous procedure, were Gaussian filtering used, would be to determine the optimal time-varying bandwidth.
The second step is to apply binary regression to the data from each trial separately to estimate Pr(t). Binary regression, within the framework of generalized linear models (McCullagh and Nelder 1990), is the equivalent of the usual regression framework, but for binary data, here spike trains transformed to a binary sequence with 1 coding for a spike at t, and 0 coding for no spike. Although BARS could also be applied here to estimate Pr(t), we found that our procedure was not sensitive to knot number and placement. We therefore obtain smooth estimates of all Pr(t) by using regression splines with the same knots used to smooth P(t). Letting r(t) denote the estimates thus obtained, then according to Eq. 1, estimates of gr(t) are 
r(t) = r(t)/(t).
Because there are often relatively few spikes per trial, there are relatively few data for fitting r(t), so that the deviations r(t) are highly variable. To reduce their variability, our next step is to identify a suitable low-dimensional representation of these deviations. This could be done very simply, for example, by smoothing further the deviations r(t), relying on splines having a small number of knots (e.g., one to three knots), or equivalently on Gaussian filtering with a large bandwidth. According to this scheme the trial-to-trial variability is allowed to be curvilinear, but is constrained to be slowly varying, across peristimulus time. The resulting smoothed deviations 
r(t) are less variable than r(t) and thus yield less variable estimates of Pr(t), as demonstrated in Fig. 1.
Our preferred method, however, is a low-dimensional principal-components representation, where the trial-to-trial fluctuation in firing-rate intensity will be assumed to have the form
 | (2) |
j(t) functions are suitably chosen curves. A special case of Eq. 2 is the constant excitability model, in which log gr(t) = w0r, discussed, for example, in Brody (1999a,b). Note that we model log gr(t) rather than gr(t) as a linear combination of the principal components to ensure that the resulting fitted gr(t) is a positive function. By restricting the form and number J of the j(t) functions we obtain an interpretable low-dimensional representation of the trial-to-trial variation, in which the coefficients wjr may be estimated from the limited data available per trial. As discussed in the following text, the j(t) functions will be taken to be principal components of the trial-to-trial variation. The set of weights wr = (w0r, w1r, w2r, ... , wJr) then describes the variability specific to trial r more simply and succinctly than would smooth functions of time.
To fit Eq. 2, we begin with the vectors of values of r(t) evaluated along the time axis: if there are time points t1, ... , tmax, then for each r we obtain the vector r(t1), ... , r(tmax). From these R vectors we compute the principal components of their covariance matrix 
. Let us write these R principal components as 1(t), ... , R(t). We then need to determine how many principal components to use, say J, and obtain the coefficients {wjr} in Eq. 2. This will effectively project each trial's deviation from the PSTH onto the space of principal components. To carry this out, a binary regression is again performed, except that gr(t) in Eq. 1 is replaced by exp[w0r + 
j=1J wjrj(t)]. Stepwise significance tests are used to determine how many of the principal components J are actually needed to represent the trial-to-trial variability. Specifically, to determine whether the contribution of some of the principal components is negligible, sequentially nested regression models should be considered: first the regression model Eq. 1 with gr(t) = 1, then the model involving trial-dependent constants with gr(t) = exp(w0r), then the model involving trial-dependent constants and the first principal component with gr(t) = exp[w0r + w1r1(t)], and so forth, up to the model involving the constants, and all principal components. These models can be fit sequentially and the deviance difference compared with a chi-squared distribution with R degrees of freedom, where R is the number of trials. Principal components may be included sequentially until the deviance difference is no longer statistically significant.
All smooth fits are obtained from generalized linear model software, available in most commercial statistical analysis packages. These regression and smoothing methods are discussed in many sources (e.g., Hastie and Tibshirani 1990; McCullagh and Nelder 1990). As we showed in Ventura et al. (2005a), spline-based generalized linear models can also be used while taking account of non-Poisson firing behavior. Such methods are applied to the simulated non-Poisson spike trains in RESULTS.
Bootstrap excursion test of time-varying dependency
Suppose two neurons labeled 1 and 2 are recorded simultaneously across multiple trials, and the task is to assess the time-varying dependency between their spike trains. Ventura et al. (2005a) defined the quantity
 | (3) |
) is the probability that neuron 1 will spike at time t and neuron 2 will spike at time t + , and P1(t) and P2(t) are the firing rates for neurons 1 and 2. Ventura et al. (2005a) then demonstrated the use of a Bootstrap procedure to assess any apparent deviation of 
(t) from 1, which would suggest extra correlation between the two neurons above the correlation induced by modulations in the firing rate.
When excitability and/or latency effects are present, we want to assess whether the neurons are correlated above what is expected from modulations in the firing rates, and from the correlation induced by trial-to-trial variability effects. To take account of trial-to-trial variability, the procedure is modified by replacing Eq. 3 with
 | (4) |
) is the probability for the rth trial that neuron 1 will spike at time t and neuron 2 will spike at time t + , and Pr1 and Pr2 are the firing rates for neurons 1 and 2 on trial r. Equation 4 assumes that the excess joint spiking activity does not itself depend on the trial. To obtain an estimate of (t), we again replace Pr1(t) and Pr2(t + ) with r1(t) and r2(t + ) as estimated earlier in Low-dimensional representation of time-varying excitability effects, and Pr12(t + ) with a smooth estimate of the joint spike count on trial r, obtained with BARS or another suitable smoothing method. The estimate of (t) is then taken to be the average of the estimates of Eq. 4 over trials r.
The Bootstrap significance test now proceeds as in Ventura et al. (2005a). In short, it is based on a computation of null bands for (t), obtained under the null hypothesis of independence. When 
(t) is above or below these bands there is potential evidence of excess or diminished joint spiking activity. If we were assessing a joint spiking activity at a single value of time t, and if (t) were above or below the null bands then there would be some evidence of excess or diminished joint activity at time t (and lag ). However, because we are examining a large number of time values, a global assessment is needed.
The only difference between the Bootstrap test of Ventura et al. (2005a) and this bootstrap test is in the sampling of Bootstrap samples. Bootstrap samples should be stochastically similar to the observed spike trains, but they should also conform to the hypothetical reality imposed by the hypothesis under investigation. Here, we are testing whether the neurons are independent above the dependency induced by modulations of firing rates and trial-to-trial variability effects. Therefore Bootstrap samples must contain excitability effects comparable to those in the observed spike trains. The modified Bootstrap procedure is as follows.
r1(t) and r2(t), where r runs through the set of trials sampled in step 1a. This is a bootstrap sample.
(t) as described above but based on this bootstrap sample rather than on the observed data, for each of interest. (t). (t).
To determine the bands accurately hL and hU, we have used bootstrap sample sizes of N = 1,000. To reduce the computation effort, a normal approximation can be used with N 50, as described in Ventura et al. (2005a).
Note that step 1 is an implementation of a parametric bootstrap. Because the trials are not exchangeable as a result of trial-to-trial variability effects, there is no nonparametric bootstrap alternative.
To perform the test of independence, we define Gobs to be the largest area of any contiguous portion of 
(t) that exceeds the bands defined above. [See Ventura et al. (2005a) for a mathematical definition of Gobs. Our Gabs is similar in spirit to k' of Ellaway and Murphy (1985).] To calculate its Bootstrap P-value, let (n)(t), n = 1, ... , N stand for the estimate of (t) obtained from the nth bootstrap sample,and Gboot(n) denote the largest area of any contiguous portion of (n)(t) that exceeds the bands defined above. We then calculate the P-value
 |
Adjustment of JPSTH and cross-correlogram for trial-to-trial variability
The joint peristimulus time histogram (JPSTH) and cross-correlogram provide useful visual representations of the correlated activity for a pair of neurons. To take account of firing-rate variation in the neurons, Aersten et al. (1989) proposed a normalized version of the JPSTH. Both the normalized JPSTH and the cross-correlogram assume stationarity across trials. In the presence of excess trial-to-trial variability, this may be violated. Once we have estimates of the trial-specific firing rates Pri(t), however, it is straightforward to adjust the JPSTH and cross-correlogram to remove the general excitability effects. The corrected JPSTH of Aertsen et al. (1989) has (t, t + ) pixel equal to
 | (5) |
(t) taken to be the value of J(t) when the spiking probabilities P1(t), P2(t), and joint spiking probability P12(t, t + ) are replaced with their observed-data counterparts, i.e., the spike counts divided by the number of trials (and the pixel width). Large values of |(t)| are evidence that the two neurons are correlated at lag .
When excitability effects are present, values of |(t)| are inflated, which may lead to the false claim that the two neurons are dependent. For values of |(t)| close to zero to still suggest independence, we define the (t, t + ) pixel of the corrected JPSTH to be
 |
(t) is the JPSTH based only on trial r, with estimate r12(t, t + ) – r1(t)r2(t + ).
A related adjustment has been used by Baker et al. (2001). The corrected cross-correlogram at lag may be then defined as
 | (6) |
, although it could be without substantially changing the methodology. Figures 4and 6 show instances of cross-correlograms adjusted for excitability effects.
|
|
In previous subsections the fitted firing rate functions r1(t) and r2(t) were not restricted to have any relationship to each other, although it is plausible that the two neurons might have either the same trial-to-trial effects or proportional trial-to-trial effects arising from common inputs. These possibilities may also be examined within the framework we have developed here, using standard generalized linear model methodology.
We will write the equality and proportionality cases as gr1(t) = gr2(t) and gr1(t) = 
· gr2(t) for all t, where is a scalar constant. A further special case is the constant excitability model gri(t) = exp(w0ri), where w0ri is a constant, for which equality and proportionality become w0r1 = w0r2 and w0r1 = log + w0r2. As explained in many texts on linear regression (the generalized linear model context being analogous), fitting these models and testing them against one another is straightforward. Such tests are performed in the next section, with results recorded in Table 4.
|
|
|
RESULTS |
|---|
|
Simulated data: Poisson spike trains
We simulate Poisson spike train data with firing rate Pr(t) given by Eq. 1 under three scenarios: no trial-to-trial variation (neuron A), constant trial-to-trial variation (neuron B), and time-varying trial-to-trial variation (neuron C). Specifics of the firing rates are summarized in Table 1 and shown in Fig. 2. All three simulated neurons have 60 trials of spike trains, and there are 200 time bins for each trial.
|
|
2 distribution with degrees of freedom equal to the difference between the degrees of freedom in the two models (60 – 0 = 60). Because the P-value is much larger than 0.05 we would conclude, correctly, that no trial-to-trial variation is present for neuron A. For neuron B the results indicate strong evidence for trial-to-trial variability (P < 10–5), but no evidence for time-dependent trial-to-trial variability (P = 0.59 > 0.05), and we would therefore correctly conclude that there is only constant trial-to-trial variability for neuron B. Figure 3 displays the true and the fitted firing rate functions for a representative subset of trials, graphically demonstrating the good fit of the constant excitability model for neuron B. The fitted rates for all 60 trials can be found in the APPENDIX, Fig. A2. Note that the no trial-to-trial variation model would have fitted the same firing rate for all trials. For neuron C the results indicate strong evidence for time-dependent trial-to-trial variability that is captured by a single principal component (P < 10–5) but no evidence that a second principal component is required to describe this variability (P > 0.05). Figure 3 displays, for a representative subset of trials, the true firing rate functions, the incorrect constant-excitability fits, and the correct time-varying fits. The fitted rates for all 60 trials can be found in the APPENDIX, Fig. A3.
|
|
Data were also simulated from a non-Poisson neuron (neuron D) that followed an IMI model (Kass and Ventura 2001), for which the conditional firing rate intensity depends not only on time t, as for Poisson processes, but also on s(t), the time of the last spike previous to time t, according to P[t, s(t)] = 
1(t)2[s(t)]. Neuron D had R = 32 trials and 300 time bins. Figure 2 shows the functions gr(t)1(t) for all 32 trials, and 2(s); specifics are in Table 1.
The estimation procedure discussed earlier in Low-dimensional representation of time-varying excitability effects was followed (as in the Poisson case), with deviances and P-values recorded in Table 2, and the resulting fitted firing rate functions for six representative trials are displayed in Fig. 3. The fits follow the true firing-rate functions on each trial quite well. The fitted rates for all 32 trials can be found in the APPENDIX, Fig. A4.
Simulated data: adjustment of dependency assessment
We simulated data from pairs of neurons in two situations. In the first scenario the two neurons (neurons E1 and E2) had independent spike trains, and both had trial-to-trial variation and latency effects; details are given in Table 3. In the second scenario the two neurons (F1 and F2) were similar to neurons E1 and E2, but they also had an excess of synchronous spikes, described by 0(t) in Table 3 and shown in Fig. 4. The firing rates of these four neurons are not shown because they are similar to neuron C in Fig. 2.
|
detected the presence of latency effects in the firing rates of the four neurons, each with values P < 10–5. We estimated the latency effects as in Ventura (2004), and shifted the trials according to the latency estimates. The methods of subsection Low-dimensional representation of time-varying excitability effects were then applied to the four neurons separately. We found that both neurons in the pair E1/E2 had firing rates suitably described by two PCAs, whereas both neurons in the pair F1/F2 had firing rates suitably described by just one PCA. To assess whether the excitability effects are shared by the two neurons, the test introduced in Trial-to-trial variation shared across neurons was used. The P-values given in Table 4 lead to the correct conclusion that the same latency and excitability effects should be fitted to each pair of simulated neurons.
Figure 4 shows the cross-correlograms and the fitted 0(t) for E1/E2 and F1/F2, with Bootstrap bands, before and after correction for latency and excitability effects.
For neurons E1/E2 (no synchrony), after the latency effect is corrected in Fig. 4B, when the excitability effect still exists, most of the cross-correlogram values are still outside the 95% bootstrap band. After both the latency and excitability effects are corrected in Fig. 4C, the cross-correlogram indicates that the two neurons are independent. Figure 4E shows the esti-mated 0(t) of the joint spiking model in Eq. 4. Before being corrected for the latency effect, the estimated 0(t) has a peak outside the Bootstrap band. The area outside the band has P-value of <0.0001. After the latency effect only is screened out, the magnitude of 0(t), shown in Fig. 4E, is reduced but with P-value of the largest area outside band equal to 0.001; all P-values are recorded in Table 5. After both the latency and excitability effect are corrected, the magnitude of 0(t), shown in Fig. 4F, is further reduced and the largest area outside the bands has a P-value of 0.63. This leads correctly to the conclusion that the two neurons are independent. The same procedure also yields the correct conclusion for the dependent pair F1/F2.
|
We assessed the effect of trial-to-trial variability on the correlation of two neurons recorded simultaneously in the primary visual cortex of an anesthetized macaque monkey (Aronov et al. 2003; units 380506.s and 380506.u, 90 deg spatial phase). The data consist of 64 trials shown in Fig. 5. The stimulus, identical for all trials, consisted of a standing sinusoidal grating that appeared at time 0, and disappeared at 237 ms.
|
5.6(t) is displayed in Fig. 6C, along with 95% Bootstrap bands; 5.6(t) exceeds the upper band for most values of t, suggesting that lagged synchrony between the two neurons is significant for almost for the entire duration of the experiment.
To adjust for trial-to-trial variability we applied the latency test of Ventura (2004) and the tests for excitability effects described above. Before doing so a preliminary check on the Poisson assumption was carried out. Based on a mean versus variance plot, there appeared to be a small amount of extra variability above that predicted by the Poisson assumption. The Poisson test of Brown et al. (2002), based on the time rescaling theorem for Poisson processes, also confirmed a very slight departure from Poisson. These effects may have been a result of trial-to-trial variation and were not large enough to warrant the more complicated procedures required for non-Poisson spike trains. When these analyses were repeated after removal of the trial-to-trial effects (described below), the Poisson assumption was judged to be adequate.
The test for latency effects in Ventura (2004) suggested that no latency effects were present, with P-values of 0.096 and 0.174 for the two neurons, respectively. However, the tests for excitability effects were highly statistically significant (P = 0.00014 and P = 0.00019, respectively). Consideration of shared excitability models showed that the most appropriate model was the shared nonconstant excitability model (P < 0.00001) against the constant excitability model (P = 0.752) for different vs. same excitability effects). Applying the sequential fitting procedure we then determined that two principal components should be used to capture the nonconstant excitability effects, with each principal component explaining 83 and 14%, respectively, of the variability in the observed pairs of spike trains. Therefore the firing rate of neuron i = 1, 2 on trial r, Pri(t), can be summarized as
 | (7) |
1(t) and 2(t) are displayed in Fig. 7, A and B.
|
 | (8) |
(t) = 5 – 1(t) – 2(t) is also plotted in Fig. 7C. That is, the firing rate for a particular trial is modulated by the function in Fig. 7C, times a gain factor.
Figure 6B shows the corrected cross-correlogram with 95% Bootstrap bands. We conclude that the correlation apparent in the uncorrected cross-correlogram (Fig. 6A) is attributed entirely to excitability effects. Figure 6C displays the estimates of 5.6(t) (corresponding to the second bin to the right of 0 in the cross-correlograms) before and after the adjustment for excitability. Before excitability adjustment, the P-value for synchrony at lag 5.6 ms is <10–5, whereas after the adjustment (t) lies within the Bootstrap bands. There is thus no evidence of synchrony at that lag, beyond the correlations induced by rate and excitability modulations. We reached similar conclusion at all lags; Fig. A5 in the APPENDIX displays the same as Fig. 6C for all lags.
Based all the analysis in this section, we conclude that 1) there is very strong evidence of time-varying excitability for these two neurons, 2) the excitability effects appear to be shared by the two neurons, and 3) there is no evidence of spike timing synchronization between these neurons.
|
|
DISCUSSION |
|---|
|
One simplification of the method used here is that we have ignored variability introduced by first smoothing P(t) and Pr(t) when we subsequently make statistical inferences based on the coefficients of the principal components. More complicated Bayesian or Bootstrap methods could provide more comprehensive inference procedures. Similarly, it would be possible to incorporate estimation of latency into a general procedure, rather than to rely on a preliminary assessment using the method of Ventura (2004).
Kass and Ventura (unpublished observations) show that excitability effects described by Eq. 1 produce spike count correlations that grow as the spike count interval increases. For realistic firing rates spike count correlations that are very small on timescales of tens of milliseconds become nonnegligible on scales of hundreds of milliseconds. This may help explain correlations that have been reported and discussed (Shadlen and Newsome 1998).
Excess trial-to-trial variability may be common in brain regions involved in higher-order processing. One way to relate trial-to-trial variation in neuronal activity to behavioral measures, such as reaction time, is to record large numbers of trials and then aggregate them into groups according to levels of the behavioral response, as in Hanes and Schall (1996). The methods used here, however, are applicable with much smaller numbers of trials.
The methods presented here have been applied to two neurons, but the general approach may, in principle, be extended to allow consideration of shared effects among many neurons recorded simultaneously.
|
|
GRANTS |
|---|
|
|
|
|
|
|
|
|
FOOTNOTES |
|---|
1 If 
(t) is an estimate of f(t), the MISE of is E {
[(t) – f(t)]2dt}, where the expectation, which is with respect to the data used to produce (t), was approximated by a sample average across many simulated samples. 
2 The choice of 2.8 ms is based on a limitation in spike recording during the experiment: if the spike of one neuron occurred immediately (less than roughly 1.3–2 ms) after the firing of the other neuron, it could not be recorded. Because it is impossible to detect a joint spike for a time lag <2.8 ms, we analyze time lags >2.8 ms or < –2.8 ms.
Address for reprint requests and other correspondence: V. Ventura, Department of Statistics and Center for the Neural Basis of Cognition, Carnegie Mellon University, Pittsburgh, PA 15213 (E-mail: vventura{at}stat.cmu.edu)
|
|
REFERENCES |
|---|
|
Aronov D, Reich DS, Mechler F, and Victor JD. Neural coding of spatial phase in V1 of the macaque monkey. J Neurophysiol 89: 3304–3327, 2003.
Azouz R and Gray CM. Cellular mechanisms contributing to response variability of cortical neurons in vivo. J Neurosci 19: 2209–2223, 1999.
Baker SN and Gerstein GL. Determination of response latency and its application to normalization of cross-correlation measures. Neural Comput 13: 1351–1377, 2001.[CrossRef][Web of Science][Medline]
Baker SN, Spinks R, Jackson A, and Lemon RN. Synchronization in monkey motor cortex during a precision grip task: I. Task-dependent modulation in single-unit synchrony. J Neurophysiol 85: 869–885, 2001.
Barbieri R, Quirk MC, Frank LM, Wilson MA, and Brown EN. Non-Poisson stimulus-response models of neural spike train activity. J Neurosci Methods 105: 25–37, 2001.[CrossRef][Web of Science][Medline]
Ben-Shaul Y, Bergman H, Ritov Y, and Abeles M. Trial to trial variability in either stimulus or action causes apparent correlation and synchrony in neuronal activity. J Neurosci Methods 111: 99–110, 2001.[CrossRef][Web of Science][Medline]
Brody CD. Correlations without synchrony. Neural Comput 11: 1537–1551, 1999a.[CrossRef][Web of Science][Medline]
Brody CD. Disambiguating different covariation. Neural Comput 11: 1527–1535, 1999b.[CrossRef][Web of Science][Medline]
Brown EN, Barbieri R, Ventura V, Kass RE, and Frank LM. A note on the time-rescaling theorem and its implications for neural data analysis. Neural Comput 14: 325–346, 2002.[CrossRef][Web of Science][Medline]
DiMatteo I, Genovese CR, and Kass RE. Bayesian curve fitting with free-knot splines. Biometrika 88: 1055–1073, 2001.
Ellaway PH and Murthy KSK. The source and distribution of short-term synchrony between gamma-motoneurones in the cat. Q J Exp Physiol Cogn Med Sci 70: 233–247, 1985.
Grün S, Riehle A, and Diesmann M. Effect of across trial non-stationarity on joint-spike events. Biol Cybern 00: 000–000, 2003.
Hanes DP and Schall JD. Neural control of voluntary movement initiation. Science 274: 427–430, 1996.
Hastie T and Tibshirani R. Generalized Additive Models. London: Chapman & Hall, 1990.
Johnson D. Point process models of single-neuron discharges. J Comput Neurosci 3: 275–299, 1996.[CrossRef][Web of Science][Medline]
Kass RE and Ventura V. A spike-train probability model. Neural Comput 13: 1713–1720, 2001.[CrossRef][Web of Science][Medline]
Kass RE, Ventura V, and Brown EN. Statistical issues in the analysis of neuronal data. J Neurophysiol 94: 8–25, 2005.
Kass RE, Ventura V, and Cai C. Statistical smoothing of neuronal data. Network: Comput Neural Syst 14: 5–15, 2005.[CrossRef]
McCullagh P and Nelder JA. Generalized Linear Models (2nd ed.). New York: Chapman & Hall, 1990.
Nawrot M, Aertsen A, and Rotter S. Single-trial estimation of neuronal firing rates. J Neurosci Methods 94: 81–91, 1999.[CrossRef][Web of Science][Medline]
Pauluis Q and Baker SN. An accurate measure of the instantaneous discharge probability, with application to unitary joint-event analysis. Neural Computation 12: 647–669, 2000.[CrossRef][Web of Science][Medline]
Shadlen MN and Newsome WT. The variable discharge of cortical neurons: implications for connectivity, computation, and information coding. J Neurosci 18: 3870–3896, 1998.
Ventura V. Testing for, and estimating latency effects for Poisson and non-Poisson spike trains. Neural Comput 16: 11: 2323–2350, 2004.[CrossRef][Web of Science][Medline]
Ventura V, Cai C, and Kass RE. Statistical assessment of time-varying dependency between two neurons. J Neurophysiol 94: 2940–2947, 2005.
Woody CD. Characterization of an adaptive filter for the analysis of variable latency neuroelectric signals. Med Biol Eng 5: 539–553, 1967.[CrossRef][Web of Science]
This article has been cited by other articles:
|
|
 |
|
  B. M. Yu, J. P. Cunningham, G. Santhanam, S. I. Ryu, K. V. Shenoy, and M. Sahani Gaussian-Process Factor Analysis for Low-Dimensional Single-Trial Analysis of Neural Population Activity J Neurophysiol, July 1, 2009; 102(1): 614 - 635. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
S. Grun Data-Driven Significance Estimation for Precise Spike Correlation J Neurophysiol, March 1, 2009; 101(3): 1126 - 1140. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 A. Ghoshal, P. Pouget, M. Popescu, and F. Ebner Early Bilateral Sensory Deprivation Blocks the Development of Coincident Discharge in Rat Barrel Cortex J. Neurosci., February 25, 2009; 29(8): 2384 - 2392. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
A. Amarasingham, T.-L. Chen, S. Geman, M. T. Harrison, and D. L. Sheinberg Spike Count Reliability and the Poisson Hypothesis J. Neurosci., January 18, 2006; 26(3): 801 - 809. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
V. Ventura, C. Cai, and R. E. Kass Statistical Assessment of Time-Varying Dependency Between Two Neurons J Neurophysiol, October 1, 2005; 94(4): 2940 - 2947. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| Visit Other APS Journals Online |
TOP
ABSTRACT
INTRODUCTION
