Signal-independent timescale analysis (SITA) and its application for neural coding during reaching and walking

Introduction

The timescales of neural encoding is important for understanding neural processing and for successfully interpreting recorded neural activity (Shadlen and Newsome, 1994; de Ruyter van Steveninck et al., 1997; Borst and Theunissen, 1999; Jacobs et al., 2009). In the context of rate-coding, we focus on characterizing those timescales, independent of which signals modulate the rate. Those may include: (i) encoded signals under experimental control, (ii) encoded task-relevant signals that are not under experimental control and thus may vary from trial to trial, and (iii) task-irrelevant signals. The last two groups are referred to as “hidden” signals, and we are especially interested in the second group of hidden task-relevant signals. Those may include, for example, the estimated state or estimation error during reaching movements (Desmurget and Grafton, 2000; Wolpert and Ghahramani, 2000; Krigolson and Holroyd, 2006; Shadmehr and Krakauer, 2008).

The potential contribution of hidden signals to trial-to-trial variability in neuronal responses has been noted in the context of different tasks, including skilled reaching movements (Churchland et al., 2006a,b; Mandelblat-Cerf et al., 2009) motor adaptation (Mandelblat-Cerf et al., 2009) and decision making (Churchland et al., 2011). During skilled reaching movements, trial-to-trial variations in preparatory neural activity were predictive of variations in reach, and were interpreted to reflect variations in movement planning (Churchland et al., 2006a). During novel visuo-motor tasks, trial-to-trial variability in spiking activity increased in early adaptation stages before returning to initial levels at the end of learning, and was interpreted to reflect enhanced exploration or adjustments of the internal models (Mandelblat-Cerf et al., 2009). During decision making, trial-to-trial variations in the underlying rate were shown to increase until the decision was made, in agreement with the hypothesis that they encode trial-specific evidence accumulation (Churchland et al., 2011).

Ideally we should quantify and analyze the contributions of the first two groups of signals, i.e., only those that are task-relevant, whether under experimental control or hidden. However, this is hampered by the inability to control or measure hidden task-relevant signals. Current methods for analyzing neuronal timescales, which have focused mainly on sensory neurons, are based on averaging the response over repeated trials with the same stimuli (Borst and Theunissen, 1999; Warzech and Egelhaaf, 1999; Butts et al., 2007), or computing the mutual information between known stimuli and their reconstruction (Rieke et al., 1997; Dimitrov and Miller, 2000). In either case, the analysis accounts only for signals that are under experimental control, i.e., the first group of signals.

Neurons in the motor system are known to encode a large number of signals (Georgopoulos, 2000; Johnson et al., 2001; Scott, 2003), and are also expected to encode hidden signals, such as the estimated state and errors (Desmurget and Grafton, 2000; Wolpert and Ghahramani, 2000; Krigolson and Holroyd, 2006; Shadmehr and Krakauer, 2008). Timescale analysis that is based on synchronized averaging would miss the contribution of those signals that are not under experimental control, including hidden task-relevant signals. This is illustrated in Figure 1A (detailed in Section Rate-modulations Signal-to-noise ratio of DSPP), where the rate is assumed to encode both a movement related signal and a hidden signal: synchronized averaging captures the rate-encoded movement but not the rate-encoded hidden signal. Instead, we focus on analyzing the timescales associated with the underlying rate modulations, independent of which signals modulate the rate. The proposed method considers asynchronous sequences of random reaching movements (Figure 1B, upper panel) so the resulting rate (Figure 1B, lower panel) is a stationary process whose variance (dashed black) captures the variance of both of the encoded signals (dashed black lines, Figure 1A upper panel and Figure 1B lower panel, respectively). Thus, in contrast with standard averaging methods, the proposed method captures also the contribution of hidden signals. Admittedly, the hidden signals may include not only task-relevant (group 2) but also task-irrelevant (group 3) hidden signals. Hence, the method is best used in comparing the timescales across different task conditions. Here we demonstrate the power of this method by comparing the timescales in neural activity during different phases of experiments with brain-machine interfaces (BMIs) and during locomotion at different speeds.

FIGURE 1

Figure 1. Illustration of the averaging approach (A) vs. the proposed approach (B) when the rate encodes both the movement and a hidden signal. Averaging methods are based on repeating the same movement (A, upper panel) and synchronizing the resulting neural activity and hence the rate (A, lower panel). Synchronized average (solid black) is considered the signal, while the variance of the hidden signal (dashed black) is considered noise. The proposed method considers asynchronous sequences of random reaching movements (B, upper panel) so the resulting rate (B, lower panel) is a stationary process whose variance (dashed black) captures the variance of both of the encoded signals [dashed black lines, (A) upper panel and (B) lower panel, respectively]. Each movement has a minimum jerk profile over 500 ms. Random reaching movements (B) are 2 s long sequences of reaching movements between random targets. Colored traces: representative samples. Black solid and dashed lines: ensemble mean and variance computed across an ensemble of 500 samples.

Timescales are assessed under the sole assumption that spike trains are realizations of dead-time modified doubly stochastic Poisson processes (DSPP) (Snyder, 1975; Cox and Isham, 1980; Johnson, 1996; Gabbiani and Koch, 1998). DSPPs are the simplest point processes that can describe rate modulations by stochastic signals. This model is pertinent for describing spike trains recorded during reaching and walking when the rate might be modulated by a number of biologically relevant stochastic signals. These signals may include not only measurable variables, such as the velocity of movement, but also internal or hidden signals (e.g., internal state estimations; Wolpert and Ghahramani, 2000; Shadmehr and Krakauer, 2008), which cannot be measured or controlled directly (Zacksenhouse et al., 2007). The effect of dead-time is explicitly investigated to account for absolute refractory period. The effects of other deviations from the DSPP assumption are demonstrated via simulations.

Within the DSPP model, the modulated rate is considered as the signal conveyed by the unit. Only rate-conditioned spike count variations are considered noise (Churchland et al., 2011). Thus, we define the rate-modulations signal-to-noise ratio (RM-SNR) as the ratio between the variance of the rate and the variance of the Poisson noise (Section Rate-modulations Signal-to-noise ratio of DSPP). Theoretical analysis, presented in Section Bandwidth Effect and Appendix A, indicates that at short bin-widths RM-SNR should increase linearly with the bin-width, while at long bin-widths it should saturate. The analysis is extended in Section Dead-time Effect and Appendix B to include the effect of refractory period (dead-time modified DSPP). Section Simulations describes the simulations that are used to demonstrate the proposed methods. Section Timescales Analysis proposes to capture the trade-off between RM-SNR and update-rate by the ratio between RM-SNR and bin-width. The variations of this ratio with the bin-width are related to the interaction between the refractory properties of the unit and the spectral properties of the encoded signals, and the bin-widths at which this ratio peaks are related to the timescales of neural encoding. The effects of deviations from the DSPP assumption are evaluated by analyzing simulated realizations of doubly stochastic Gamma processes (DSGPs).

After demonstrating the analysis on simulated realizations of both DSPPs and doubly stochastic Gamma processes (DSGPs), it is applied to investigate the timescales of the neural activity recorded from cortical motor units during two experiments with Monkeys (described briefly in Section Experimental Methods). Sections SNR During BMI Experiments and Timescale During BMI Experiments present the analysis of spike-trains recorded during experiments with BMIs and demonstrate that the average timescale agrees well with the experimentally selected bin-width, though the timescales of individual units vary. The effect of the refractory period is assessed in Section Refractory Effects. Most importantly, Section BMI and Training Effects demonstrates that the timescales did not change significantly when switching to brain control or with training. During locomotion, the analysis identified the units whose timescales varied with the speed of walking—in agreement with the spectral analysis (Section Timescales During Walking). We conclude in Section Conclusions and Discussion, by discussing the proposed signal-independent timescale analysis (SITA) and its significance for investigating the timescales of neural-rate coding under different task conditions, and how they are related to the spectral content of the encoded signals and the refractory properties of the unit.

Materials and Methods

Modeling and Analysis Methods

Rate-modulations Signal-to-noise ratio of DSPP

The Poisson process is the simplest point process in which the probability of an event—a spike in the case of neural activity, is assumed to be independent of the history of the spike train. Inhomogeneous Poisson processes are the simplest point processes that can describe rate modulations, i.e., the probability of a spike may change with time (Snyder, 1975; Cox and Isham, 1980; Johnson, 1996; Gabbiani and Koch, 1998). Here we consider DSPP, which is a subclass of inhomogeneous Poisson processes whose rate is modulated by stochastic signals (Snyder, 1975; Cox and Isham, 1980; Johnson, 1996; Gabbiani and Koch, 1998). Thus, the statistics of spike trains generated by DSPPs are determined by two factors: stochastic changes in instantaneous spike-rate and Poisson probability of spike occurrence.

Considering long intervals of movements, the stochastic signals that modulate the instantaneous rate, and thus the instantaneous rate itself, are assumed to be stationary. This is justified since we consider long sequences of movements (reaching or stepping) that are not synchronized to any event, and thus may start at an arbitrary phase of the movement. This approach is illustrated in Figure 1 and contrasted with synchronized averaging, using a simple example where the rate is assumed to encode both the position (along 1-dimension) during a reaching movement and an independent hidden signal. Figure 1A illustrates the averaging method, where the same reaching movement (upper panel) is repeated, and the resulting spike-rates (lower panel, depicting a sample of 5 realizations) are synchronized to movement initiation. The resulting synchronized rate is a non-stationary process with time-varying mean (solid black line, averaged over a simulated sample of 500 realizations). Figure 1B illustrates the proposed approach, showing a sample of five movement sequences (upper panel) starting at random phases (for illustration, each sequence includes four consecutive reaching movements between random targets). The resulting spike-rates (lower panel) are not synchronized and thus generate a stationary process with constant mean and variance (solid and dashed black lines, respectively, based on a sample of 500 realizations). The rate-variance (Figure 1B, lower panel) captures the contributions of both the movement (whose variance is depicted as black dashed line in Figure 1B, upper panel), and the hidden signal (whose variance is depicted as black dashed line in Figure 1A, lower panel).

Thus, the instantaneous spike-rate λ(t) = λ₀ +$\tilde{\lambda }$ (t) is assumed to be a stationary stochastic process with mean rate λ₀ and zero-mean rate modulations $\tilde{\lambda }$ (t). In realizations of DSPP, the distribution of spike-counts, N_T, in bins of size T is determined by the integrated spike-rate during the bin:

$\begin{array}{cc}{\Lambda }_T(t)={\displaystyle {\int }_{t-T/2}^{t+T/2}\lambda }(\sigma )d\sigma ,& (1)\end{array}$

and its statistics are given by (Snyder, 1975; Zacksenhouse et al., 2007):

$\begin{array}{ll}\begin{array}{l}E[N_T]=E[{\Lambda }_T]\\ Var[N_T]=Var[{\Lambda }_T]+E[{\Lambda }_T]\end{array}\hfill & (2)\hfill \end{array}$

The last equation can be interpreted as a decomposition of the total variance in the spike-counts into the variance of the modulated rate, Var[Λ_T], and the variance of the noise E[Λ_T ] = E[N_T] (Zacksenhouse et al., 2007). A similar decomposition was derived for general point processes from the total variance equation (Churchland et al., 2011), however the restriction to DSPPs assures that the meaning of the two sources of the neural variance is clearly defined (as further detailed in Section Conclusions and Discussion), and facilitates explicit analysis of the effect of the refractory period, as detailed in Section Dead-time Effect.

Considering the modulated rate as the signal, the RM-SNR is defined as:

$\begin{array}{ll}SN{R}_T\equiv \frac{Var[{\Lambda }_T]}{E[{\Lambda }_T]},\hfill & (3)\hfill \end{array}$

Using Equation (2), RM-SNR can be estimated from the mean and variance of the binned spike-counts as:

$\begin{array}{ll}S\widehat{N}{R}_{N_T}=\frac{Var[N_T]-E[N_T]}{E[N_T]}=F-1\hfill & (4)\hfill \end{array}$

where the hat denotes estimation, and F is the Fano factor, defined as the ratio between the variance and the mean of the spike-counts (Dayan and Abbott, 2001). Thus, for the simple DSPP case, the RM-SNR is the deviation of the Fano-factor from its value for the homogeneous Poisson process, where F = 1. Positive RM-SNR reflects rate modulations and characterizes irregular spike trains (i.e., spike trains whose variance is larger than their mean and thus larger than the variance of a homogeneous Poisson process with the same mean). Expressing the SNR is advantageous for further analysis of the timescales, as detailed in Section Timescales Analysis.

Bandwidth effect

The relationship between RM-SNR and the spectrum S_{$\tilde{\lambda }$} (ω) of the instantaneous rate $\tilde{\lambda }$ (t) is analyzed in Appendix A. Equation (A1) shows that:

$\begin{array}{ll}SN{R}_T=\frac{1}{{\lambda }_{0}T}\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }{\left|G_{W_T}(\omega )\right|}^2}{S}_{\tilde{\lambda }}(\omega )d\omega \hfill & (5)\hfill \end{array}$

where $\left|G_{W_T}(\omega )\right|=T\left(\frac{\sin (\omega T/2)}{\omega T/2}\right)$ is the Fourier transform of the rectangular window of duration T (Bendat and Peirsol, 2000). Furthermore, under the assumption that the spectrum of the instantaneous spike-rate S_{$\tilde{\lambda }$} (ω) is band limited in the range [−ω_max, ω_max], Appendix A derives two asymptotes for RM-SNR that hold at short and long bin-widths, respectively. For short bin-widths T, i.e., when ω_max «2π /T, Equation (A2) indicates that RM-SNR increases linearly with the bin-width. For long bin-widths T, i.e., when ω_max »2π /T, Equation (A3) indicates that RM-SNR saturates.

Dead-time effect

During absolute refractory period, or dead-time τ_d, the instantaneous firing rate is zero regardless of the modulating signals. The effect of the dead-time on the estimated RM-SNR is detailed in Appendix B, and, assuming that the dead-time is short compared to the dynamics of the instantaneous rate, an exact expression for S$\widehat{N}$R_NT is derived (Equation B4).

To gain more insight into the effect of the dead-time, the estimated S$\widehat{N}$R_NT is approximated (Appendix B, Equation B9) by:

$\begin{array}{ll}S\widehat{N}{R}_{N_T}({\tau }_d)\cong \frac{SN{R}_T}{\left[1+{\tau }_d{\lambda }_0\right]}-\frac{{\tau }_d{\lambda }_0\left(2+{\tau }_d{\lambda }_0\right)}{{\left[1+{\tau }_d{\lambda }_0\right]}^2}\hfill & (6)\hfill \end{array}$

The approximation indicates that S$\widehat{N}$R_NT estimated from a realization of a dead-time modified DSPP is a scaled and downward shifted version of SNR_T associated with the original, dead-time free, rate. The scaling reflects the effect of the dead-time on the integrated rate (Equation B6), and the shift reflects the additional effect on the variance of the spike-counts (Equation B3).

Simulations

The analysis is demonstrated on simulated realizations of point processes with known characteristics. Specifically, the instantaneous rate was generated as a stationary stochastic process by: (a) passing a unit variance white Gaussian noise through a 2nd order Butterworth filter with known cut-off frequency f_cut, (b) scaling the resulting signal to achieve the desired variance of the instantaneous rate, and (c) adding a constant mean rate λ₀. Finally, simulated spike trains were generated from the resulting instantaneous rate using the inverse distribution function technique (Johnson, 1996). Spike trains were generated as realizations of DSPPs, doubly stochastic Gamma processes (DSGP; Fujiwara et al., 2009) with constant shape parameter κ, or dead-time modified DSPPs or DSGPs with fixed dead-time. In summary, the simulated spike trains were characterized by five parameters: bandwidth defined by the cut-off frequency f_cut of the low-pass filter, mean (λ₀) and variance (V_λ) of the instantaneous rate, dead-time τ_d, and the shape parameter for DSGP (κ, κ = 1 for DSPP).

We note that for homogeneous Gamma processes with κ >1, the probability of firing after a spike is reduced, resulting in relative refractory period, and the variance of the spike-counts is lower than for the homogeneous Poisson process with the same mean. In contrast, for homogeneous Gamma processes with κ <1, the probability of firing after a spike is enhanced, and the variance of the spike counts is larger than for the homogeneous Poisson process with the same mean.

Timescales analysis

While SNR_T, and hence S$\widehat{N}$R_NT, increase with the bin-width, the update rate decreases (Wu et al., 2006) and the reaction time increases (Cohen and Newsome, 2009). We quantify the trade-off between S$\widehat{N}$R_NT and update-rate by evaluating their product, i.e., the ratio S$\widehat{N}$R_NT/T, as a function of the bin-width T, and characterize the timescale of neural encoding as the bin-width at which this ratio saturates or peaks. This is demonstrated in Figure 2, which depicts the variations of S$\widehat{N}$R_NT/T with the bin-width for simulated realizations of DSPPs and DSGPs generated from the same realization of the instantaneous rate. The instantaneous rate was generated with mean λ₀ = 15 s⁻¹, variance V_λ = 48 s⁻², and cut-off frequency f_cut = 1.0 Hz. The realizations analyzed in Figure 2A are dead-time free, while the realizations analyzed in Figure 2B include dead-time τ_d = 1 ms.

FIGURE 2

Figure 2. Timescale analysis of simulated realizations of doubly stochastic Poisson and Gamma processes (DSPP, DSGP) (A) and dead-time modified (τ_d = 1) DSPP and DSGP (B). DSGPs are simulated with κ = 0.85, κ = 1 (i.e., DSPP), and κ = 1.15. All simulations were generated from the same instantaneous rate function generated with mean λ₀ = 15 s⁻¹, variance V_λ = 48 s⁻², and cutoff frequency f_cut = 1.0 Hz. Estimated S$\widehat{N}$R_NT/T (Equation 4) are compared with the values computed from the rate [using Equation (3) for the dead-time free DSPP in (A), and Equation (B4) for the dead-time modified DSPP in (B)]. Analysis is conducted at bin-widths of 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 175, 200, 250, 300, 500, 750, and 1000 ms.

Timescale analysis of DSPPs. The estimated S$\widehat{N}$R_NT/T from both the dead-time free DSPP (Figure 2A, blue line) and dead-time modified DSPP (Figure 2B, blue line) closely match the values computed from the rate (Equations 3 and B4, respectively, black lines). The curves depict the expected effects of the spectral content and dead-time. At long bin-widths, S$\widehat{N}$R_T should reach a constant asymptote (Equation A3), so SNR_T/T should indeed decrease with the bin-width. The additive bias in estimating S$\widehat{N}$R_NT from dead-time modified DSPP is constant (Equation 6), independent of the bin-width, so, even when the spike-trains include a refractory period, the effect on the estimated S$\widehat{N}$R_NT/T would diminish at long bin-widths, as is evident in Figure 2B. At the other extreme, i.e., when the bin-widths are short enough for the linear asymptote (Equation A2) to hold, SNR_T/T should be approximately constant, in agreement with the estimated S$\widehat{N}$R_NT/T for DSPP (blue curve, Figure 2A). Hence, the short bin-widths at which S$\widehat{N}$R_NT/T saturates are indicative of the timescales above which the slower than linear increase in SNR does not justify the reduction in update rate.

Dead-time introduces a constant bias in estimating S$\widehat{N}$R_NT (Equation 6), which becomes increasingly significant as SNR_T diminishes at short bin-widths. Since the bias is negative, the estimated S$\widehat{N}$R_NT/T becomes smaller as the bin-width becomes shorter, in agreement with the short bin-width trend of the blue curve in Figure 2B. The effect of the dead-time at short bin-widths limits the benefit of reducing the bin-width, and results in a peak in the S$\widehat{N}$R_NT/T curve. The bin-width at which this peak occurs maximizes the trade-off between S$\widehat{N}$R_NT and update rate. We also note that below this bin-width the distortion due to the dead-time increases significantly.

Timescale analysis of DSGPs. Applied to doubly stochastic Gamma processes (DSGPs) with scale parameter (κ) larger/smaller than 1, the estimated S$\widehat{N}$R_NT/T under/over estimates the value computed from the rate, as demonstrated in Figure 2. In realizations of Gamma processes with κ >1, the occurrence of a spike has an inhibitory effect, resulting in relative refractory period. The estimated S$\widehat{N}$R_NT/T curve in such a case (red line, Figure 2A) or when an absolute refractory period is also included (red line, Figure 2B) depict a peak similar to that for dead-time modified DSPP. As noted above (Section Simulations) realizations of homogeneous Gamma processes with κ >1 are characterized by spike-count variance that is smaller than the mean. The positive estimated S$\widehat{N}$R_NT (which implies variance larger than mean, see Equation 4) reflects the contribution of the variance of the stochastic rate.

In realizations of DSGPs with κ <1, the occurrence of a spike enhances, rather than inhibits, the probability of firing. The estimated S$\widehat{N}$R_NT/T curve in such a case (green line, Figure 2A) increases sharply as the bin-width become shorter without saturating. For the specific parameters selected here, the excitatory recovery after a spike dominates the inhibitory effect due to the dead-time, and a similar increase is observed even with dead-time (Figure 2B).

Timescales from SNR_NT/T curves. In summary, the analysis and simulations indicate that: (i) DSPPs are characterized by decreasing S$\widehat{N}$R_NT/T curves that saturate at short bin-widths. In those cases, the bin-widths at which the S$\widehat{N}$R_NT/T curve saturates are indicative of the timescales above which the slower than linear increase in SNR does not justify the reduction in update rate. (ii) dead-time modified DSPP or DSGP with relative refractory period (κ >1) are characterized by S$\widehat{N}$R_NT/T curves that peak at some bin-width. The bin-width at which the S$\widehat{N}$R_NT/T peaks maximizes the trade-off between S$\widehat{N}$R_NT and update rate, while limiting the distortion due to the refractory effect. (iii) DSPGs with excitatory recovery (κ <1) are characterized by decreasing S$\widehat{N}$R_NT/T curves that remain convex even at short bin-widths without saturating. In those cases, the appropriate time-scale is not evident from the proposed analysis, since the higher S$\widehat{N}$R_NT/T at short bin-widths reflect recovery effects rather than actual rate modulations. Further analysis is needed to select the bin-width that optimizes the trade-off between high S$\widehat{N}$R_NT/T and small distortion due to the excitatory recovery effect.

Variance, bandwidth and dead-time effects. The effects of the variance, bandwidth and dead-time are further evaluated in Figure 3, based on the approximation derived in Equation 6. Figure 3A establishes that the approximated S$\widehat{N}$R_NT/T (dashed black line) captures well the main effect of the dead-time. Hence, the above effects are evaluated by approximating S$\widehat{N}$R_NT directly from the simulated stochastic rate functions. Unless otherwise specified, the nominal parameters of the stochastic rate functions are: mean rate λ₀ = 15 s⁻¹, variance V_λ = 48 s⁻², bandwidth f_cut = 1 Hz, and dead-time τ_d = 0.7 ms.

FIGURE 3

Figure 3. Approximated S$\widehat{N}$R_NT/T for dead-time modified DSPP (Equation 6) in comparison with the S$\widehat{N}$R_NT/T computed from the rate (Equation B4) or estimated from the spike-counts (A). The approximation is used to demonstrate the effects of the variance (B) bandwidth (C) and dead-time (D). Unless otherwise specified, the nominal parameters of the stochastic rate functions are: mean rate λ₀ = 15 s⁻¹, variance V_λ = 48 s⁻², bandwidth f_cut = 1 Hz, and dead-time τ_d = 0.7 ms.

Figure 3B demonstrates that the major effect of increasing the variance (from V_λ = 24 s⁻² to V_λ = 48 s⁻² and V_λ = 72 s⁻²) is to increase the SNR. This is also accompanied by a shift in the peak location toward shorter bin-widths, due to the interaction between SNR_T and dead-time (Equation 6). Figure 3C demonstrates that as the cutoff frequency f_cut increases from 0.5 Hz to 1.0 Hz, and 1.5 Hz, S$\widehat{N}$R_NT/T peaks at shorter bin-widths, capturing well the reduction in the relevant timescales. Figure 3D demonstrates that as the dead-time increases to 0.7 and 1.4 ms, S$\widehat{N}$R_NT/T peaks at longer bin-widths. This trend reflects the increasing effect of the dead-time on the magnitude of the negative bias in Equation (6) (see note following Equation B9). Thus, for dead-time modified DSPP, the bin-width at which S$\widehat{N}$R_NT/T peaks reflects the conflicting effects of the bandwidth of the encoded signals and the refractory period.

Experimental Methods

The analysis tools developed above to estimate the RM-SNR of spike-trains and evaluate their timescales are applied to spike-trains recorded form cortical units in two experiments with Monkeys: (a) BMI experiments and (b) Locomotion experiments.

Ethics statement

All experiments were conducted with approved protocols from the Duke University Institutional Animal Care and Use Committee and were in accordance with the NRC/NIH guidelines for the Care and Use of Laboratory Animals.

BMI experiments

As detailed in Carmena et al. (2003), the BMI experiments were conducted with macaque monkeys and included three control modes: (i) pole control, (ii) brain control with hand movements (BCWH), and (iii) brain control without hand movements (BCWOH). During pole control the monkey controlled the cursor on the screen using a hand held pole. Data from pole control was used to train a linear filter to predict the velocity from the neural activity. During brain control the cursor was controlled to follow the velocity predicted from the recorded neural activity using the previously trained linear filter. Initially, the monkey continued to move its hand during brain control (BCWH), but starting at the 7th session, the monkey stopped moving the hand (BCWOH). Spike trains were recorded from single and multi-units (see Section Refractory Effects for more details) in the primary motor area (M1), the pre-motor area (PMd), the supplementary motor area (SMA), and somato-sensory area (S1). Here we focus on significantly modulated units (n = 163), i.e., those units for which the null hypothesis that their ISI distribution was generated from a homogeneous Poisson process can be rejected at α = 0.05 significance level (Zacksenhouse et al., 2007).

Locomotion experiments

During the locomotion experiments, rhesus macaques either stood or walked bipedally on a treadmill in three different speeds, 12.5, 25, and 50 cm/s, as detailed in Fitzsimmons et al. (2009). Spike trains were recorded from single and multi-units (66 and 34%, respectively) in regions associated with the representation of the lower limbs in the primary motor (M1) and somatosensory (S1) cortices.

Results

SNR during BMI Experiments

Figures 4A,C depict the mean S$\widehat{N}$R_NT and mean normalized S$\widehat{N}$R_NT estimated from spike-trains of significantly modulated cortical units (n = 163) recorded during a BMI experiment, as a function of the bin-width T. Assuming ergodicity, Equation (4) was used to estimate the S$\widehat{N}$R_NT of each unit from its spike-counts statistics in non-overlapping windows of 2 min. The mean S$\widehat{N}$R_NT in Figure 4A was derived by first computing ensemble mean (over n = 163 units) at each window, and then computing the mean and standard deviation over 10 consecutive windows in each of the experimental modes (pole control and brain control with and without hand movements). Thus, error bars in Figure 4A depict the standard deviations of the ensemble-means over the 10 windows of time. Variations across the ensemble of units are described in Figure 4C, which depicts the ensemble mean and standard deviations of the normalized S$\widehat{N}$R_NT. The normalized S$\widehat{N}$R_NT for each unit was computed by first averaging across all windows, and then normalizing by a scaling factor. The scaling factor was selected so the normalized S$\widehat{N}$R_NT at bin-width of 300 ms is the same for all units and equals the original ensemble mean. Thus, error bars in Figure 4C depict the standard deviations of the normalized S$\widehat{N}$R_NT across the ensemble of significantly modulated units (resulting in zero standard deviation at 300 ms, which was used for normalizations).

FIGURE 4

Figure 4. Mean estimated S$\widehat{N}$R_NT (A) and normalized S$\widehat{N}$R_NT (C), and the corresponding S$\widehat{N}$R_NT/T (B,D) of neural activity recorded from significantly modulated units (n = 163) during BMI experiments. BMI experiments included pole control, brain control with hand movements (WH), and without hand movements (WOH). S$\widehat{N}$R_NT was computed from the spike-train of each unit in non-overlapping windows of 2-min. (A,B) are based on first averaging over the units at each window, and then computing the mean and standard deviation over 10 consecutive windows in each of the experimental mode. Hence, error-bars in (A,B) mark the standard deviations across 10 consecutive windows. (C,D) are based on first averaging the S$\widehat{N}$R_NT of each unit across all windows, and then scaling it so the normalized S$\widehat{N}$R_NT at bin-width of 300 ms is the same for all units and equal the original ensemble mean. Ensemble mean and standard deviation of the normalized S$\widehat{N}$R_NT were computed across units. Hence, error bars in (C,D) mark the standard deviations across the ensemble of units (resulting in zero standard deviation at 300 ms, which was used for normalizations). Dashed lines in (A) depict the linear fits for short bin-widths, over the range 50–150 ms.

Figures 4A,C demonstrate that the mean estimated S$\widehat{N}$R_NT increases approximately linearly at short bin-widths and approaches a constant level at long bin-widths, in agreement with the expected asymptotes (Equations A2 and A3). Another striking feature is that at each bin-width, the mean S$\widehat{N}$R_NT in brain control is always higher than that in pole control; and is highest in brain control without hand movements. This extends our previous results (Zacksenhouse et al., 2007) for T = 100 ms—the bin-width used for operating the BMI.

Timescale during BMI Experiments

Figures 4B,D demonstrate that the mean S$\widehat{N}$R_NT/T and mean normalized S$\widehat{N}$R_NT /T increase at short bin-widths, reach a wide peak around 100 ms, and decrease at long bin-widths. This pattern agrees well with the analysis and approximation of the S$\widehat{N}$R_NT/T of dead-time modified DSPP, and with estimations derived from the simulated spike-trains of dead-time modified DSPPs (or DSGPs with κ <1, with or without dead-time) depicted in Figure 2.

Focusing on individual units from PMd and M1, the S$\widehat{N}$R_NT /T curves shown in Figure 5 demonstrate two typical patterns: (i) curves with a wide peak around 100 ms (Figures 5C,D), similar to the average curves in Figures 4B,D and to the curves for dead-time modified DSPP or DSGP with κ <1, in Figure 2, or (ii) decreasing curves (Figures 5A,B), similar to the curves for DSPP or DSGP with κ >1 in Figure 2. The second type could be further divided depending on whether the curve saturates at short bin-widths as for DSPPs, or remains convex as for DSGPs with κ >1.

FIGURE 5

Figure 5. Estimated S$\widehat{N}$R_NT /T as a function of the bin-width for four representative cortical units in PMd (A,C) and M1 (B,D), during pole control, and brain control with and without hand movements, WH and WOH, respectively.

Units were classified as exhibiting a wide peak around 100 ms if the peak S$\widehat{N}$R_NT/T was in the range of 50–150 ms, and as exhibiting decreasing S$\widehat{N}$R_NT/T if the peak was below 50 ms. Note that S$\widehat{N}$R_NT/T of few units peaked at bin-widths longer than 150 ms, and were not included in either group. To assure proper peak-detection, we restricted this classification and further analysis of the timescales to the significantly modulated units whose maximum S$\widehat{N}$R_NT/T during pole control was larger than 0.5 s⁻¹ (n = 109). Of those units, 59.7, 64.2, and 61.5% exhibited a wide peak around 100 ms during pole control, BCWH and BCWOH, respectively, while 25.7, 21.1, and 22.0% exhibited decreasing curves, respectively. The second group was further divided depending on whether the decreasing curve was convex at short bin-widths, i.e., whether u₂ < (u₁ + u₃)/2 where u₁, u₂, u₃ are the mean S$\widehat{N}$R_NT/T in the bin-width intervals 30–50, 60–80, and 90–100 ms, respectively. About 15% of the curves were convex at short bin-widths (17.4, 11.0, and 15.6%, in pole control BCWH and BCWOH, respectively). In those cases, the time-scale analysis is not conclusive, and is either not appropriate since the recovery function is dominated by excitatory effects (as clarified in Section Timescales Analysis) or should be extended to shorter bin-widths.

The S$\widehat{N}$R_NT/T estimated form spike-trains recorded during BMI experiments (Figures 4B,D, 5) vary with the bin-width in a similar way to the S$\widehat{N}$R_NT/T of simulated spike-trains (Figure 2). Decreasing curves, which saturate at short bin-widths (Figure 5B, pole and BCWH) agree well with the assumption that the spike-trains are realizations of DSPP while decreasing curves that remain convex at short bin-widths (Figure 5A, pole and BCWH) agree well with the assumption that the spike-trains are realizations of DSGP with κ <1. Curves that peak (Figures 5C,D) agree well with the assumption that the spike trains are realizations of dead-time modified DSPP or DSGP with κ >1 (with or without dead-time), which are characterized by relative, rather than absolute, refractory period.

Refractory effects

To assess the potential contribution of refractory effects, we estimated the refractory interval from the ISI distribution. The minimum ISI in all units was 0.675 ms, and most units (95% of the significantly modulated units) had minimum ISI of less than 1 ms. Since a single ISI is prone to measurement errors, we estimated the dead-time as the interval which was exceeded by 99% of the ISIs in at least one mode of operation. Of the 109 units described above, 50% had dead-time greater than 1.6 ms, and thus could be considered single-units (Hatsopoulos et al., 2004). The percent of units having decreasing curves was similar, independent of whether the dead-time was short or long (below or above 1.6 ms); but the percent of units having time-scales longer than 150 ms was larger for units with long dead-time (22.4% on average, across all modes, compared with 8% for units with short refractory interval). Furthermore, the bin-widths at which S$\widehat{N}$R_NT /T peaked had a moderate positive correlation with the product of the dead-time and the mean rate (coefficient of correlation of 0.4, 0.42, 0.48 in pole control BCWH and BCWOH, respectively). Both the positive correlation and its moderate value are in agreement with Equation (6), which indicates that this product reduces the S$\widehat{N}$R_NT (and hence shifts the peak of the S$\widehat{N}$R_NT /T toward longer bin-widths, as clarified in Section Timescales Analysis), but that S$\widehat{N}$R_NT is also affected by other variables (i.e., SNR_T).

BMI and Training Effects

Figure 5 indicates that the bin-width at which the S$\widehat{N}$R_NT/T peaks remains approximately the same even after switching to brain control. Specifically, the S$\widehat{N}$R_NT/T curves in pole control, BCWH and BCWOH peak at: 30, 40, and 30 ms for the PMd unit in panel (A), 130, 150, and 140 ms for the PMd unit in panel (B); 30, 40, and 30 ms, for the M1 unit in panel (C); and 150, 150, and 110 ms for the M1 unit in panel (D). The distribution of the bin-widths at which the estimated S$\widehat{N}$R_NT/T peak and how they change when switching to brain control is summarized in Figure 6. The bin-width at which the estimated S$\widehat{N}$R_NT/T peaks during pole control is similar to (within 50 ms of) the bin-width at which it peaks during brain control for most of the units (80.7 and 67.9% for BCWH and BCWOH, respectively).

FIGURE 6

Figure 6. Distribution of timescales during BMI experiments. Comparison of bin-widths at which S$\widehat{N}$R_NT/T curves peak during pole control and during either brain control with hand movements, (WH, A) or brain control without hand movements (WOH, B). Only significantly modulated units whose peak S$\widehat{N}$R_NT/T during pole control was larger than 0.5 s⁻¹ were included (n = 109).

The spectral content of the spike-counts (in 100 ms bins) of the four units whose timescales were analyzed in Figure 5, are depicted in Figure 7. Specifically, the power spectrum density (PSD) was computed from the spike-counts in bins of 100 ms after subtracting the mean and expressed in dB. The main effect of the transition to brain control is the increase in power, which is higher in brain control than in pole control and is highest in brain control without hand movements. This is consistent with the higher RM-SNR in brain control compared to pole control. The overall bandwidth of the activity of each unit remains relatively similar across the different BMI control modes [even when, as in panel (D), the PSD of the activity during pole control portrays a small peak which is surpassed by the higher overall power during brain control]—consistent with the relatively similar shapes and peak locations of the different S$\widehat{N}$R_NT/T curves for each unit in Figure 5.

FIGURE 7

Figure 7. Power spectrum density (PSD) of the neural activity analyzed in Figure 4. The spectrums were computed from spike-counts in 100 ms bins, with frequency resolution of 0.039 Hz.

Figure 8 demonstrates that the shape of the S$\widehat{N}$R_NT/T curves is also invariant to training. At later sessions, as training progresses, the SNR at each control mode and bin-width decreases (as detailed in Zacksenhouse et al., 2007) for bin-width of 100 ms). However, the effect of bin-width is similar and the S$\widehat{N}$R_NT/T curves peak at a similar bin-width in all sessions.

FIGURE 8

Figure 8. Mean estimated S$\widehat{N}$R_NT/T across the population of significantly modulated cortical units recorded in three different BMI sessions during pole control (PC), and brain control with and without hand movements (BCWH and BCWOH, repectively) as a function of the bin-width. The three sessions include: an early session (3rd day); a mid session (7th day), and a late session (10th day). The 7th day is the first one in which BCWOH was performed, and in subsequent sessions BCWH was not perfomrned.

Timescales during Walking

Figures 9, 10 depict the estimated S$\widehat{N}$R_NT/T and the spectrum of representative cortical units in the motor leg area while the monkey was standing (black) or walking at increasing speeds (12.5, 25, and 50 cm/s). The S$\widehat{N}$R_NT/T curves in Figures 9A,B peak at progressively shorter bin-widths as the speed of walking increases. During standing, the peak S$\widehat{N}$R_NT /T is shallower and appears at longer bin-widths than during walking. Figures 10A,B indicate that the spectrums of the spike-counts (in 100 ms bins) of these spike-trains are dominated by single spectral lines that appear at higher frequencies as the speed of walking increases. During standing, the dominating spectral line appears at very low frequencies and is smaller compared to the spectral lines during walking. Thus, this simple case, in which the spectrums are dominated by single spectral lines, helps to demonstrate the expected effect of changes in the spectral content of the modulating signals on the timescale. However, the timescale emerges from the interaction between the spectral content and the refractory properties of the individual unit. Hence, even though the dominating spectral lines appear in the same frequencies (Figures 10A,B) the timescales are unit specific (Figures 9A,B).

FIGURE 9

Figure 9. Estimated S$\widehat{N}$R_NT/T as a function of the bin-width for four representative cortical units (A–D) in the motor leg area while the monkey is standing (black) or walking at increasing speeds (12.5, 25, and 50 cm/s).

FIGURE 10

Figure 10. Power Spectrum Density (PSD) of the neural activity analyzed in Figure 9. PSDs were computed from the binned neural activity with 100 ms bins and frequency resolution of 0.039 Hz.

In contrast, the neural activity from other units in the leg area are characterized by S$\widehat{N}$R_NT/T curves that either decrease without peaking (Figure 9D) or depict a peak that does not shift as expected with the speed of walking (Figure 9C). In both cases the neural activity during standing is characterized by higher S$\widehat{N}$R_NT/T than during walking. The corresponding spectrums (Figures 10C,D) indicate that the power of the spike-counts recorded from these two units is indeed larger during standing than during walking. The spectral lines associated with the speed of walking (at the frequencies in which they appear in the spectrums depicted in panels A and B) are very small. Finally, for both units the neural activity recorded during running has more power in the low-frequency range than during walking. This may explain the unexpected shift in the peak of the S$\widehat{N}$R_T /T toward longer rather than shorter bin-widths when the Monkey is running compared to walking (Figure 9C). The convex shape of the S$\widehat{N}$R_NT/T curve in Figure 9D even at short bin-widths suggests that the activity of this unit is dominated by excitatory recovery effects, and the timescale is not conclusive. Nevertheless, since the S$\widehat{N}$R_NT/T curves are similar during walking and running (independent of the speed) but much higher during standing, it can be concluded that signals associated with the frequency of walking have little effect on the spike-rate of these two units, and that their activity is more strongly modulated during standing.

Additional analysis (not shown) indicates that the mean S$\widehat{N}$R_NT/T and mean normalized S$\widehat{N}$R_NT/T across significantly modulated units with peak above 0.5 s⁻¹ (n = 48) depict a wide peak around 50–110 ms, in most cases. The only exception is the normalized S$\widehat{N}$R_NT/T during standing, which is highest at the shortest bin-width analyzed (25 ms). Around the peak, the magnitude of the S$\widehat{N}$R_NT/T are similar independent of the speed so the curves overlap each other. At longer bin-widths S$\widehat{N}$R_NT/T during slow walking is higher than during fast walking or running, and is highest during standing.

Conclusions and Discussion

SNR Estimation

We suggest a computational framework for estimating the SNR in spike trains without relying on any assumption about the encoded signals. The suggested approach captures all the encoded signals, including hidden and internal signals, and not only those controlled experimentally. The proposed approach has two limitations: (i) it accounts also for possible task irrelevant rate modulations, and (ii) the estimation method is derived under the assumption that the spike-trains are realization of DSPP. The first limitation is a necessary side-effect of capturing the hidden or internal signals. The second limitation is partially relaxed by considering dead-time modified DSPP and doubly stochastic Gamma processes.

The dependency of RM-SNR on bin-width was first analyzed for DSPP and two asymptotes were derived for bin-widths that are either very short or very long compared with the inverse of the bandwidth. In particular, it was shown that RM-SNR increases linearity at short bin-widths and reaches saturation at long bin-widths. Extending the analysis to dead-time modified DSPP, it was shown that the estimated RM-SNR is reduced by both scaling and shifting factors that do not depend on the bin-width. The bias effect becomes more significant as the RM-SNR becomes smaller at shorter bin-widths. Using simulations of doubly-stochastic Gamma processes, we also evaluate other recovery effects in which the occurrence of a spike temporarily and partially inhibits (as in relative refractory period) or enhances the probability of firing.

Inhomogeneous vs. Doubly Stochastic Point Processes

Doubly stochastic point processes are inhomogeneous point processes with stochastic rate function. The alternative model for describing rate modulations is the inhomogeneous point process with a deterministic rate function (Cox and Isham, 1980; Kass and Ventura, 2001; Nawrot et al., 2008). The deterministic model is based on the assumption that the rate function depends solely on the experimentally controlled conditions, and thus that it can be estimated from the post-stimulus time histogram. This assumption might be valid for peripheral neurons that are affected almost exclusively by the experimentally controlled stimulus, but less for cortical neurons which may encode additional hidden signals, not directly under experimental control. These signals may vary from trial-to-trial and thus can only be modeled as stochastic. During movements, motor planning (Churchland et al., 2006a,b) and prediction errors (Wolpert and Ghahramani, 2000; Shadmehr and Krakauer, 2008), for example, may vary from trial-to-trial, while in decision making the accumulated evidence (while potentially sharing a common trend), may vary from trial-to-trial (Churchland et al., 2011). In these cases, doubly stochastic point processes are essential for capturing the variance of all the signals that might be encoded in the neural activity.

Variance Decomposition

Under the assumption that spike-trains are realizations of DSPPs, the variance of the spike-counts is decomposed into two terms associated with the signal and noise. A similar decomposition was developed in Churchland et al. (2011) where the resulting two terms were described as the variance of the underlying “intensity command” and the variance of the neural activity given that command. That decomposition was applied to uncover how trial-to-trial neural variability changes during decision making and reveal the underlying nature of neural processing. While the application and some of the analysis details differ (as discussed below), we share the most important hypotheses: (i) the neural activity may represent hidden signals, i.e., signals that are not directly controlled and hence may vary from trial to trial, (ii) decomposition of the variance of neural activity facilitates estimating the variance of the underlying rate, including modulations by those hidden signals, and (iii) variations in the variance of the underlying rate with time (or task) may reveal important information about the neural processing.

The decomposition developed in Churchland et al. (2011) is based on the total variance equation, which holds for any doubly stochastic point process. The total variance is decomposed into the variance of the conditional expectation (the variance of the expected spike-count given the rate parameter and the point process variance (the expectation of the conditional variance). Under the restriction to DSPP made here, the first term equals the variance of the rate parameter, which is assumed to encode the modulating signals. Otherwise, the relationship of the first term to the variance of the rate parameter is more complex. In particular, for dead-time modified DSPPs and for DSGPs, the first term is proportional to the variance of the rate-parameter but not equal to it (see Equations B1 and B6 for the dead-time modified DSPP). Hence, the restriction to DSPP made here assures that the meaning of the first term is clearly defined. Furthermore, it enables explicit analysis of the effect of the refractory period on both terms. Another major difference is that the analysis in Churchland et al. (2011) was applied to estimate the synchronized trial-to-trial variance in the underlying rate-process (intensity command) at each time along the trial. In contrast, we applied the analysis to estimate asynchronous variance in the underlying rate process.

Timescales of Neural Activity

As bin-width increases, SNR increases but update-rate is compromised (Wu et al., 2006). This trade-off is captured by the ratio of SNR to bin-width, SNR_T/T. Thus, the bin-width at which SNR_T/T peaks optimizes the trade-off between SNR and update rate and is suggested as a criterion for characterizing the timescales of neural rate-coding.

Our analysis suggests that the estimated S$\widehat{N}$R_NT/T is affected by both the spectral content of the encoded signals and the recovery properties of the unit. The effects of the bandwidth and absolute refractory period were analyzed theoretically and demonstrated via simulations of DSPPs and dead-time modified DSPPs. Their conflicting effects give rise to a peak in the S$\widehat{N}$R_NT/T curve, which characterizes the timescale of neural rate-coding.

Deviations from the DSPP Assumption

The DSPP assumption is violated when the firing rate depends on the history of the spike train. In these cases, the estimated S$\widehat{N}$R_NT may differ from the actual SNR_T, due to history-dependent modulations of the instantaneous rate. We restrict the analysis to the simple, yet common case, of renewal point processes, in which the dependency on the history is captured by the time elapsed since the last spike. The effect of absolute refractory period was analyzed theoretically, and was shown to result in a peak in the S$\widehat{N}$R_NT/T curve. Thus, the dead-time limits the benefit of reducing the bin-width, since the estimated S$\widehat{N}$R_NT becomes increasingly distorted. The bin-width at which the peak occurs maximizes the trade-off between S$\widehat{N}$R_NT and update rate, and characterizes the timescale of neural rate-coding.

Using simulations we also demonstrated that relative refractory periods have similar effect, resulting in S$\widehat{N}$R_NT/T curves that are also characterized by a peak. So the timescales are similarly characterized by the bin-widths at which the S$\widehat{N}$R_NT/T curves peak. In contrast, when the recovery period is characterized by enhanced, rather than inhibited, probability of firing, the S$\widehat{N}$R_NT/T over-estimates (rather than under-estimates) the SNR_T/T. Thus, as the bin-width decreases, S$\widehat{N}$R_NT/T increases without peaking or saturating. In those cases, the appropriate time-scale is not conclusive from the proposed analysis, since the higher S$\widehat{N}$R_NT/T at short bin-widths reflect recovery effects rather than actual signal modulations. Further analysis is needed to select the bin-width that optimizes the trade-off between high S$\widehat{N}$R_NT/T and small distortion due to the excitatory recovery effect.

Timescales during BMI Experiments

Estimating RM-SNR from neural activity recorded during BMI experiments, we observe that while its magnitude depends on the BMI mode and on training, the timescale revealed by the peak of the S$\widehat{N}$R_NT/T curve remains the same.

The observation that RM-SNR increases when switching to brain control indicates that the variance of the rate increases—and thus that the variance of some of the encoded signals increases. The observation that the timescales remain invariant suggests that those signals are task-relevant, and we are currently investigating the hypothesis that they include state estimation and control signals, within the framework of optimal control. However, their exact nature and decoding remains an open issue.

Most of the analyzed S$\widehat{N}$R_NT/T curves (62% on average across the different control modes, n = 109) were characterized by a wide peak at 50–150 ms. For this class of curves (and for the 15% of curves that peak at longer bin-widths), decreasing the bin-width below the peak may compromise decoding, unless the effect of the refractory period is accounted for. Indeed BMI applications based on spike-counts used bin-widths in the range between 30 and 100 ms (Wessberg et al., 2000; Serruya et al., 2002; Kim et al., 2008; Velliste et al., 2008), consistent with the range of timescales revealed here. In particular, the BMI experiments analyzed here were conducted with 100 ms bins (Carmena et al., 2003) within the wide peak of the estimated S$\widehat{N}$R_NT/T curve depicted in Figure 4B. Our analysis suggests that this selection provides a good trade-off between the SNR and update rate.

Curves whose peak was below 50 ms were divided into those that were convex or concave at short bin-widths (15 and 8% on average, respectively). Concave shape at short bin-widths is expected from S$\widehat{N}$R_NT/T curves derived from DSPPs, where this ratio should saturate at short bin-widths as detailed in Section Timescales Analysis (and Figure 2A). In this case, the range of bin-widths at which the curve saturates is indicative of the relevant timescales for neural decoding. At longer bin-widths, the slower than linear increase in SNR does not justify the reduction in update rate. Convex S$\widehat{N}$R_NT/T are expected when the unit has excitatory recovery period (in which the probability of firing is enhanced). As detailed above, the timescale analysis in this case is not conclusive.

Since the relevant timescales vary across individual units, it might be beneficial to apply different bin-widths when binning different units. This may explain the potential advantage of multi-resolution binning (Kim et al., 2005). Recent decoding methods are based on point process filters that operate at 5 ms bins (Shanechi et al., 2013a,b). The BMI tested in Shanechi et al. (2013a) was based on multi-unit recording, where the dead-time effect might be negligible, and hence the timescales could be short. Furthermore, point process filters are based on estimating the conditional probability of the spike-counts given the rate, which is computationally more efficient when the bin is small (and the number of spikes that may occur with significant probability is limited). Hence, the choice of small bins may reflect computational constraints. The analysis conducted here suggests that unless recovery effects are negligible or accounted for, very short bin-widths may compromise the trade-off between SNR and update-rate.

Timescales during Locomotion Experiments

The overall shape of S$\widehat{N}$R_NT/T curves estimated from spike trains recorded during walking is similar to those derived from the spike trains recorded during the BMI experiments. However, the peak location in the locomotion experiments varied with the task, at least for some units, and shifted to lower bin-widths when the speed of walking increased. This is in agreement with the expected change in the bandwidth of encoded signals associated with the experimentally controlled speed of walking. Indeed, units whose timescales vary in this way are characterized by PSDs dominated by a single peak at the frequency of walking. Thus, the locomotion test-case demonstrates the expected effect of changes in the spectral content of the modulating signals on the timescales. However, timescales depend not only on the spectral content but also on the refractory properties of the units. So even if the spectral content is similar the timescales are unit-specific.

In summary, we developed a method for assessing the SNR and timescales of neuronal activity independent of the encoded signals and observe that many units depict well-defined timescales of the order of tens to hundreds msec. The timescale emerges from the interaction of the recovery properties of the unit and the spectral content of the modulating signals. In particular, the variations of the estimated S$\widehat{N}$R_NT/T at short bin-widths is related to the recovery properties of the unit. During a BMI experiment, the average timescale across the population agrees well with the experimentally selected bin-width, though individual units depict different timescales. Furthermore, while SNR increases after switching to brain control and decreases with training, timescales remain similar. During locomotion, the analysis identified units whose timescales varied consistently with the experimentally controlled speed of walking, though the specific timescale reflected also the recovery properties of the unit. Hence, the proposed method provides a hypothesis free tool for investigating possible contributions of hidden signals, not under experimental control, to neural modulations and the effect of task conditions on their magnitude and timescales.

Conflict of Interest Statement

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

Acknowledgments

This research was supported by the Ashkenazy Handicap Research Fund and TASP 2018172 to Miriam Zacksenhouse and by the Anne W. Deane Endowed Chair Fund, Telemedicine and Advanced Technology Research Center (TATRC) W81XWH-08-2-0119 to Miguel A. L. Nicolelis, DARPA, the James S. McDonnel Foundation, NIH and NSF to Miguel A. L. Nicolelis.

References

Bair, W., Koch, C., Newsome, W., and Britten, K. (1994). Power spectrum analysis of bursting cells in area MT in the behaving monkey. J. Neurosci. 74, 2870–2892.

Pubmed Abstract | Pubmed Full Text

Bendat, J. S., and Peirsol, A. G. (2000). Random Data Analysis and Measurements Procedures. New York, NY; London; Sydney, NSW; Toronto, ON: John Wiley and Sons, Inc.

Borst, A., and Theunissen, F. E. (1999). Information theory and neural coding. Nat. Rev. 2, 947–957.

Pubmed Abstract | Pubmed Full Text

Butts, D. A., Weng, C., Jianzhing, J., Yeh, C. I., Lesica, N. A., Alonso, J. M., et al. (2007). Temporal precision in the neural code and the timescales of natural vision. Nature 446, 92–96. doi: 10.1038/nature06105

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Carmena, J. M., Lebedev, M. A., Crist, R. E., O'Doherty, J., and Santucci, D. M. (2003). Learning to control a brain–machine interface for reaching and grasping by primates. PLoS Biol. 1:E42. doi: 10.1371/journal.pbio.0000042

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Churchland, A. K., Kiani, R., Chaudhuri, R., Wang, X. J., Pouget, A., and Shadlen, M. N. (2011). Variance as a signature of neural computations during decision making. Neuron 69, 818–831. doi: 10.1016/j.neuron.2010.12.037

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Churchland, M. M., Afshar, A., and Shenoy, K. V. (2006a). A central source of movement variability. Neuron 52, 1085–1096. doi: 10.1016/j.neuron.2006.10.034

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Churchland, M. M., Yu, B. M., Ryu, S. I., Santhanam, G., and Shenoy, K. V. (2006b). Neural variability in premotor cortex provides a signature of motor preparation. J. Neurosci. 26, 3697–3712. doi: 10.1523/JNEUROSCI.3762-05.2006

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cohen, M. R., and Newsome, W. T. (2009). Estimated of the contribution of single neurons to perception depend on timescales and noise correlations. J. Neurosci. 29, 6635–6648. doi: 10.1523/JNEUROSCI.5179-08.2009

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cox, D. R., and Isham, V. (1980). Point Processes. Vol. 12. Boca Raton, FL: CRC Press.

Dayan, P., and Abbott, L. F. (2001). Theoretical Neuroscience. Cambridge, MA: MIT Press.

de Ruyter van Steveninck, R. R., Lewen, G. D., Strong, S. P., Koberie, R., and Bialek, W. (1997). Reproducibility and variability in neural spike trains. Science 275, 1805–1808. doi: 10.1126/science.275.5307.1805

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Desmurget, M., and Grafton, S. (2000). Forward modeling allows feedback control for fast reaching movements. Trends Cogn. Sci. 4, 423–431. doi: 10.1016/S1364-6613(00)01537-0

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Dimitrov, A. G., and Miller, J. P. (2000). Natural time scales for neural encoding. Neurocomputing 32, 1027–1034. doi: 10.1016/S0925-2312(00)00275-7

CrossRef Full Text

Fitzsimmons, N. A., Lebedev, M. A., Peikon, I. D., and Nicolelis, M. A. L. (2009). Extracting kinematic parameters for monkey bipedal walking from cortical neuronal ensemble activity. Front. Integr. Neurosci. 3:3. doi: 10.3389/neuro.07.003.2009

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Fujiwara, K., Aihara, K., and Suzuki, H. (2009). “Interspike interval statistics obtained from non-homogeneous gamma spike generator,” in Neural Information Processing, eds C.-S. Leung, M. Lee, and J. H. Chan (Berlin; Heidelberg: Springer), 143–150. doi: 10.1007/978-3-642-10677-4_16

CrossRef Full Text

Gabbiani, F., and Koch, C. (1998). “Principles of spike train analysis methods,” in Neuronal Modeling: from Synapses to Networks, 2nd Edn., eds C. Koch and I. Segev (Cambridge, MA: MIT Press), 313–360.

Georgopoulos, A. P. (2000). Neural aspects of cognitive motor control. Curr. Opin. Neurobiol. 10, 238–241. doi: 10.1016/S0959-4388(00)00072-6

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hatsopoulos, N., Joshi, J., and O'Leary, J. G. (2004). Decoding continuous and discrete motor behaviors using motor and premotor cortical ensembles. J. Neurophysiol. 92, 1165–1174. doi: 10.1152/jn.01245.2003

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Jacobs, A. L., Fridman, G., Douglas, R. M., Alam, N. M., Latham, P. E., Prusky, G. T., et al. (2009). Ruling out and ruling in neural codes. Proc. Natl. Acad. Sci. U.S.A. 106, 5936–5941. doi: 10.1073/pnas.0900573106

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Johnson, D. H. (1996). Point process models of single neuron discharge. J. Comput. Neurosci. 3, 275–299. doi: 10.1007/BF00161089

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Johnson, M. T. V., Mason, C. R., and Ebner, T. J. (2001). Central processes for the multiparametric control of arm movements in primates. Curr. Opin. Neurobiol. 11, 684–688. doi: 10.1016/S0959-4388(01)00269-0

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kass, R. E., and Ventura, V. (2001). A spike-train probability model. Neural Comput. 13, 1713–1720. doi: 10.1162/08997660152469314

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kim, S. P., Carmena, J. M., Nicolelis, M. A. L., and Principe, J. C. (2005). “Multiresolution representation and data mining of neural spikes for brain-machine interfaces,” in IEEE Neuroeng. Conf. (Arlington, TX).

Kim, S. P., Simeral, J. D., Hochberg, L. R., Donoghue, J. P., and Black, M. J. (2008). Neural control of computer cursor velocity by decoding motor cortical spiking activity in humans with tetraplegia. J. Neural Eng. 5, 455. doi: 10.1088/1741-2560/5/4/010

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Krigolson, O. E., and Holroyd, C. B. (2006). Evidence for hierarchical error processing in the human brain. Neuroscience 137, 13–17. doi: 10.1016/j.neuroscience.2005.10.064

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Mandelblat-Cerf, Y., Paz, R., and Vaadia, E. (2009). Trial-to-trial variability of single cells in motor cortices is dynamically modified during visuomotor adaptation. J. Neurosci. 29, 15053–15062. doi: 10.1523/JNEUROSCI.3011-09.2009

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Nawrot, M. P., Boucsein, C., Rodriguez Molina, V., Riehle, A., Aertsen, A., and Rotter, S. (2008). Measurement of variability dynamics in cortical spike trains. J. Neurosci. Methods 169, 374–390. doi: 10.1016/j.jneumeth.2007.10.013

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Rieke, F., Warland, D., de Ruyter van Steveninck, R. R., and Bialek, W. (1997). Spikes: Exploring the Neural Code. Cambridge: MIT Press.

Scott, S. H. (2003). The role of primary motor cortex in goal directed movements: insights from neurphysiological studies on non-human primates. Curr. Opin. Neurobiol. 13, 671–677. doi: 10.1016/j.conb.2003.10.012

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Serruya, M. D., Hatsopoulos, N. G., Paninski, L., Fellows, M. R., and Donoghue, J. P. (2002). Instant neural control of a movement signal. Nature. 416, 141–142. doi: 10.1038/416141a

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Shadlen, M. N., and Newsome, W. T. (1994). Noise, neural codes and cortical organization. Curr. Opin. Neurobiol. 4, 569–579. doi: 10.1016/0959-4388(94)90059-0

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Shadmehr, R., and Krakauer, J. W. (2008). A computational neuroanatomy for motor control. Exp. Brain Res. 185, 359–381. doi: 10.1007/s00221-008-1280-5

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Shanechi, M. M., Orsborn, A., Gowda, S., and Carmena, J. M. (2013a). “Proficient BMI control enabled by closed-loop adaptation of an optimal feedback-controlled point process decoder,” in Transl and Comp Motor Control Meeting (San Diego, CA).

Shanechi, M. M., Williams, Z. M., Wornell, G. W., Hu, R. C., Powers, M., and Brown, E. N. (2013b). A real-time brain-machine interface combining motor target and trajectory intent using an optimal feedback control design. PLoS ONE 8:e59049. doi: 10.1371/journal.pone.0059049

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Snyder, S. L. (1975). Random Point Processes. New York, NY; London; Sydney, NSW; Toronto, ON: John Wiley & Sons, Inc.

Vannucci, G., and Teich, M. C. (1981). Dead-time modified photo-count mean and variance for chaotic radiation. J. Opt. Soc. Am. 71, 164–170. doi: 10.1364/JOSA.71.000164

CrossRef Full Text

Velliste, M., Perel, S., Spalding, M. C., Whitford, A. S., and Schwartz, A. B. (2008). Cortical control of a prosthetic arm for self-feeding. Nature 453, 1098–1101. doi: 10.1038/nature06996

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Warzech, A. K., and Egelhaaf, M. (1999). Variability in spike trains during constant and dynamic stimulation. Science 283, 1927–1030. doi: 10.1126/science.283.5409.1927

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wessberg, J., Stambaugh, C. R., Kralik, J. D., Beck, P. D., Laubach, M., Chapin, J. K., et al. (2000). Real-time prediction of hand trajectory by ensembles of cortical neurons in primates. Nature 408, 361–365. doi: 10.1038/35042582

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wolpert, M. W., and Ghahramani, Z. (2000). Computational principles of movement neuroscience. Nat. Neurosci. 3, 1212–1217. doi: 10.1038/81497

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wu, W., Gao, Y., Bienenstock, E., Donoghue, J. P., and Black, M. (2006). Bayesian population decoding of motor cortical activity using a Kalman filter. Neural Comput. 18, 80–118. doi: 10.1162/089976606774841585

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Zacksenhouse, M., Lebedev, M. A., Carmena, J. M., O'Doherty, J. E., Henriquez, C. S., and Nicolelis, M. A. L. (2007). Cortical modulations increase in early sessions with brain-machine interface. PLoS ONE 2:e619. doi: 10.1371/journal.pone.0000619

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Appendix A: Bandwidth Effect

The binning process is equivalent to passing the instantaneous spike-rate through a finite time integrator followed by down-sampling and spike-counting (see Bair et al., 1994 for a similar formalization of the recording of spike-trains with finite time resolution). Specifically, the integrated spike-rate Λ_T(t) = λ₀T + $\tilde{\Lambda }$_T (t), where ${\tilde{\Lambda }}_T(t)={\displaystyle {\int }_{t-T/2}^{t+T/2}\tilde{\lambda }}(\sigma )d\sigma$, is obtained by convolving the time-dependent component of the instantaneous spike-rate $\tilde{\lambda }$ (σ) with a rectangular window of duration T: $\tilde{\Lambda }$_T(t) = $\tilde{\lambda }$ (t) ⊗ W_T (t) [⊗ denotes convolution, and W_T (t) is a rectangle signal of duration T (Bair et al., 1994)]. The effect of down-sampling will be considered later. For now we note that the spectrum of $\tilde{\Lambda }$_T(t), denoted by S_{$\tilde{\Lambda }$} (ω), is related to the spectrum S_{$\tilde{\lambda }$} (ω) of $\tilde{\lambda }$ (t) by: S_{$\tilde{\Lambda }$} (ω) = |G_WT (ω)|² S_{$\tilde{\lambda }$} (ω), where $\left|G_{W_T}(\omega )\right|=T\left(\frac{\sin (\omega T/2)}{\omega T/2}\right)$ is the Fourier transform of the rectangular window (Bendat and Peirsol, 2000). The variance of Λ_T(t) can be computed from the spectrum S_{$\tilde{\Lambda }$} (ω) as: $Var[{\Lambda }_T(t)]=\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }{S}_{\tilde{\Lambda }}}(\omega )d\omega$, while its mean is E[Λ_T(t)] = E[N_T(t)] = λ₀T. Using Equation (3), the SNR of the binned spike-counts is given by:

$\begin{array}{ll}SN{R}_T=\frac{1}{{\lambda }_{0}T}\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }{\left|G_{W_T}(\omega )\right|}^2}{S}_{\tilde{\lambda }}(\omega )d\omega \hfill & (\text{A}1)\hfill \end{array}$

Assuming that the spectrum of the instantaneous spike-rate S_{$\tilde{\lambda }$} (ω) is band limited in the range [−ω_max, ω_max], the SNR reaches two asymptotes for short and long bin-widths. For short bin-widths T, i.e., when ω_max «2π /T, the transfer function of the rectangular window |G_WT (ω)| in the range [−ω_max, ω_max] can be approximated by |G_WT (0)| = T, so, from Equation (A1):

$\begin{array}{ll}SN{R}_T\cong \frac{1}{{\lambda }_{0}T}\frac{T^2}{2\pi }{\displaystyle {\int }_{-{\omega }_{\max }}^{{\omega }_{\max }}{S}_{\tilde{\lambda }}}(\omega )d\omega =\frac{P_{\tilde{\lambda }}}{{\lambda }_0}T,\hfill & (\text{A}2)\hfill \end{array}$

where $P_{\tilde{\lambda }}=\frac{1}{2\pi }{\displaystyle {\int }_{-{\omega }_{\max }}^{{\omega }_{\max }}{S}_{\tilde{\lambda }}}(\omega )d\omega$ is the power of the time-dependent component of the instantaneous rate. Hence, at short bin-widths, the SNR of spike-counts generated from DSPP with band-limited instantaneous rate should increase linearly with the bin-width T.

For long bin-widths T, i.e., ω_max »2π /T, the Fourier transform of the rectangular window is concentrated around zero, so only low frequencies pass the averaging process and affect the average rate. Denoting by S₀ the spectrum of the instantaneous rate at low frequencies, Equation (A1) implies that at long bin-widths the SNR saturates at:

$\begin{array}{ll}SN{R}_T\underset{T\to \infty }{\to }\frac{1}{{\lambda }_{0}T}\frac{S_0}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }{\left|G_{W_T}(\omega )\right|}^2}d\omega =\frac{S_{0}T}{{\lambda }_{0}T}=\frac{S_0}{{\lambda }_0}.\hfill & (\text{A}3)\hfill \end{array}$

Binning is usually performed in non-overlapping windows, so the binned spike-counts are related to the samples of the integrated spike-rate taken at sampling interval T. At short bin-widths, the Nyquist frequency ω_Ny ≡ ω_sampling /2 ≡ π /T satisfies the Nyquist criterion ω_Ny >ω_max. Hence, there is no aliasing, and Equation (A2) holds, implying that the SNR of binned spike counts should increase linearly with the bin-width. At long bin-widths, the Nyquist criterion is not satisfied, and the side bands of the spectrum of the rectangular window |G_WT (ω)|² fold back into the main lobe. Nevertheless, as long as the spectrum is constant over most of the side-lobes, the total power in the main lobe (including all the folded back side-lobes) is the same as the integral computed in Equation (A3) (Bair et al., 1994). Hence the SNR of binned spike counts should saturates at long bin-widths.

Appendix B: SNR of Dead-Time Modified DSPP

Assuming the dead-time is short compared to the dynamics of the instantaneous rate, the statistics of the spike-counts depend on the dead-time modified integrated rate: ${\overline{\Lambda }}_T={\displaystyle {\int }_{t-T/2}^{t+T/2}\frac{\lambda (\sigma )}{1+{\tau }_d\lambda (\sigma )}}d\sigma$, and are given by [see Equations 3 and 9 in Vannucci and Teich (1981)]:

$\begin{array}{cc}E[N_T]=E\left[{\overline{\Lambda }}_T\right]& (\text{B}1)\end{array}$

$\begin{array}{ll}E[N_T^2]=E\left[{\left({\overline{\Lambda }}_T\right)}^2\right]+E\left[{\displaystyle {\int }_{t-T/2}^{t+T/2}\frac{\lambda (\sigma )}{{\left[1+{\tau }_d\lambda (\sigma )\right]}^3}}d\sigma \right]\hfill & (\text{B}2)\hfill \end{array}$

So:

$\begin{array}{ll}Var[N_T]=Var[{\overline{\Lambda }}_T]+E\left[{\displaystyle {\int }_{t-T/2}^{t+T/2}\frac{\lambda (\sigma )}{{\left[1+{\tau }_d\lambda (\sigma )\right]}^3}}d\sigma \right]\hfill & (\text{B}3)\hfill \end{array}$

Substituting Equations (B1) and (B3) in Equation (4), the estimated SNR from a dead-time modified DSPP is:

$\begin{array}{ll}S\widehat{N}{R}_{N_T}({\tau }_d)=\frac{Var\left[{\overline{\Lambda }}_T\right]+E\left[{\displaystyle {\int }_{t-T/2}^{t+T/2}\frac{\lambda (\sigma )}{{\left[1+{\tau }_d\lambda (\sigma )\right]}^3}}d\sigma \right]-E\left[{\overline{\Lambda }}_T\right]}{E\left[{\overline{\Lambda }}_T\right]}\hfill & (\text{B}4)\hfill \end{array}$

The SNR associated with the dead-time modified integrated rate, SNR_T, can be defined in a similar way to the definition of SNR_T for the original integrated rate (Equation 3):

$\begin{array}{cc}S\overline{N}{R}_T\equiv \frac{Var\left[{\overline{\Lambda }}_T\right]}{E\left[{\overline{\Lambda }}_T\right]}& (\text{B}5)\end{array}$

For τ_d$\tilde{\lambda }$ (σ) « 1, the integral defining Λ_T can be approximated by: ${\overline{\Lambda }}_T\cong \frac{1}{1+{\tau }_d{\lambda }_0}{\displaystyle {\int }_{t-T/2}^{t+T/2}\lambda }(\sigma )d\sigma =\frac{{\Lambda }_T}{1+{\tau }_d{\lambda }_0}$, so its statistics can be approximated by:

$\begin{array}{cc}E[{\overline{\Lambda }}_T]\cong \frac{E\left[{\Lambda }_T\right]}{1+{\tau }_d{\lambda }_0}& (\text{B}6)\end{array}$

$\begin{array}{ll}Var[{\overline{\Lambda }}_T]\cong \frac{Var\left[{\Lambda }_T\right]}{{\left[1+{\tau }_d{\lambda }_0\right]}^2}\hfill & (\text{B}7)\hfill \end{array}$

Hence, the SNR_T can be related to SNR_T by:

$\begin{array}{ll}S\overline{N}{R}_T({\tau }_d)\cong \frac{SN{R}_T}{\left[1+{\tau }_d{\lambda }_0\right]}\hfill & (\text{B}8)\hfill \end{array}$

Using a similar approximation, the second term in Equations (B2–B4) can be approximated as:

$\begin{array}{ll}E\left[{\displaystyle {\int }_{t-T/2}^{t+T/2}\frac{\lambda (\sigma )}{{\left[1+{\tau }_d\lambda (\sigma )\right]}^3}}d\sigma \right]\cong \frac{E\left[{\Lambda }_T\right]}{{\left[1+{\tau }_d{\lambda }_0\right]}^3}\hfill & (\text{B}9)\hfill \end{array}$

Substituting these approximations in Equation (B4) results in an approximated relationship between the estimated SNR from a dead-time modified DSPP and the SNR of a dead-time free DSPP having the same instantaneous rate:

$\begin{array}{ll}\begin{array}{l}S\widehat{N}{R}_{N_T}({\tau }_d)\cong S\overline{N}{R}_T-\frac{{\tau }_d{\lambda }_0\left(2+{\tau }_d{\lambda }_0\right)}{{\left[1+{\tau }_d{\lambda }_0\right]}^2}\cong \frac{SN{R}_T}{\left[1+{\tau }_d{\lambda }_0\right]}\\ -\frac{{\tau }_d{\lambda }_0\left(2+{\tau }_d{\lambda }_0\right)}{{\left[1+{\tau }_{d_0}\right]}^2}\end{array}\hfill & (\text{B}10)\hfill \end{array}$

Note that the second term increases with the dead-time both absolutely and as a fraction of the first term.