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

What are the relevant timescales of neural encoding in the brain? This question is commonly investigated with respect to well-defined stimuli or actions. However, neurons often encode multiple signals, including hidden or internal, which are not experimentally controlled, and thus excluded from such analysis. Here we consider all rate modulations as the signal, and define the rate-modulations signal-to-noise ratio (RM-SNR) as the ratio between the variance of the rate and the variance of the neuronal noise. As the bin-width increases, RM-SNR increases while the update rate decreases. This tradeoff is captured by the ratio of RM-SNR to bin-width, and its variations with the bin-width reveal the timescales of neural activity. Theoretical analysis and simulations elucidate how the interactions between the recovery properties of the unit and the spectral content of the encoded signals shape this ratio and determine the timescales of neural coding. The resulting signal-independent timescale analysis (SITA) is applied to investigate timescales of neural activity recorded from the motor cortex of monkeys during: (i) reaching experiments with Brain-Machine Interface (BMI), and (ii) locomotion experiments at different speeds. Interestingly, the timescales during BMI experiments did not change significantly with the control mode or training. During locomotion, the analysis identified units whose timescale varied consistently with the experimentally controlled speed of walking, though the specific timescale reflected also the recovery properties of the unit. Thus, the proposed method, SITA, characterizes the timescales of neural encoding and how they are affected by the motor task, while accounting for all rate modulations.


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 taskirrelevant (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.
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 deadtime 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 Signalto-noise ratio of DSPP). Theoretical analysis, presented in Section Bandwidth Effect and Appendix A, indicates that at short binwidths 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 updaterate 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 signalindependent 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.

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) = λ 0 +λ(t) is assumed to be a stationary stochastic process with mean rate λ 0 and zero-mean rate modulationsλ(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: and its statistics are given by (Snyder, 1975;Zacksenhouse et al., 2007): 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: Using Equation (2), RM-SNR can be estimated from the mean and variance of the binned spike-counts as: 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λ(ω) of the instantaneous rateλ(t) is analyzed in Appendix A. Equation (A1) shows that: where G W T (ω) = T sin (ωT/2) ωT/2 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λ(ω) 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 binwidths 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 SNR N T is derived (Equation B4).
To gain more insight into the effect of the dead-time, the estimated SNR N T is approximated (Appendix B, Equation B9) by: The approximation indicates that SNR N T 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 λ 0 . 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 (λ 0 ) 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 SNR N T , 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 SNR N T and update-rate by evaluating their product, i.e., the ratio SNR N T /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 SNR N T /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 λ 0 = 15 s −1 , variance V λ = 48 s −2 , 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  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.
reach a constant asymptote (Equation A3), so SNR T /T should indeed decrease with the bin-width. The additive bias in estimating SNR N T 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 SNR N T /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 SNR N T /T for DSPP (blue curve, Figure 2A). Hence, the short bin-widths at which SNR N T /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 SNR N T (Equation 6), which becomes increasingly significant as SNR T diminishes at short bin-widths. Since the bias is negative, the estimated SNR N T /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 binwidths limits the benefit of reducing the bin-width, and results in a peak in the SNR N T /T curve. The bin-width at which this peak occurs maximizes the trade-off between SNR N T 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 SNR N T /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 SNR N T /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 SNR N T (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 SNR N T /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 N T /T curves.
In summary, the analysis and simulations indicate that: (i) DSPPs are characterized by decreasing SNR N T /T curves that saturate at short bin-widths. In those cases, the bin-widths at which the SNR N T /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 SNR N T /T curves that peak at some bin-width. The bin-width at which the SNR N T /T peaks maximizes the trade-off between SNR N T and update rate, while limiting the distortion due to the refractory effect. (iii) DSPGs with excitatory recovery (κ < 1) are characterized by decreasing SNR N T /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 SNR N T /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 SNR N T /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 SNR N T /T (dashed black line) captures well the main effect of the dead-time. Hence, the above effects are evaluated by approximating SNR N T directly from the simulated stochastic rate functions. Unless otherwise specified, the nominal parameters of the stochastic rate functions are: mean rate λ 0 = 15 s −1 , variance V λ = 48 s −2 , 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 −2 to V λ = 48 s −2 and V λ = 72 s −2 ) 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, SNR N T /T peaks at shorter binwidths, capturing well the reduction in the relevant timescales. Figure 3D demonstrates that as the dead-time increases to 0.7 and 1.4 ms, SNR N T /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 SNR N T /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 spiketrains 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 multiunits (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.

SNR DURING BMI EXPERIMENTS
Figures 4A,C depict the mean SNR N T and mean normalized SNR N T 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 SNR N T of each unit from its spike-counts statistics in non-overlapping windows of 2 min. The mean SNR N T 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 SNR N T . The normalized SNR N T 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 SNR N T 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 SNR N T across the ensemble of significantly modulated units (resulting in zero standard deviation at 300 ms, which was used for normalizations).
Figures 4A,C demonstrate that the mean estimated SNR N T 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 SNR N T 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.

Figures 4B,D demonstrate that the mean SNR N T /T and mean
normalized SNR N T /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 SNR N T /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 SNR N T /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.
Units were classified as exhibiting a wide peak around 100 ms if the peak SNR N T /T was in the range of 50-150 ms, and as exhibiting decreasing SNR N T /T if the peak was below 50 ms. Note that SNR N T /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 SNR N T /T during pole control was larger than 0.5 s −1 (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 2 < (u 1 + u 3 ) /2 where u 1 , u 2 , u 3 are the mean SNR N T /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 SNR N T /T estimated form spike-trains recorded during BMI experiments (Figures 4B,D, 5) vary with the bin-width in a similar way to the SNR N T /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 spiketrains 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. M1 (B,D), during pole control, and brain control with and without hand movements, WH and WOH, respectively.

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 timescales 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 SNR N T /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 SNR N T (and hence shifts the peak of the SNR N T /T toward longer bin-widths, as clarified in Section Timescales Analysis), but that SNR N T is also affected by other variables (i.e., SNR T ). . The distribution of the bin-widths at which the estimated SNR N T /T peak and how they change when switching to brain control is summarized in Figure 6. The bin-width at which the estimated SNR N T /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). 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 SNR N T /T curves for each unit in Figure 5. Figure 8 demonstrates that the shape of the SNR N T /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 SNR N T /T curves peak at a similar bin-width in all sessions.

TIMESCALES DURING WALKING
Figures 9, 10 depict the estimated SNR N T /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 SNR N T /T curves in Figures 9A,B peak at

FIGURE 6 | Distribution of timescales during BMI experiments.
Comparison of bin-widths at which SNR N T /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 SNR N T /T during pole control was larger than 0.5 s −1 were included (n = 109).

Frontiers in Computational Neuroscience
www.frontiersin.org August 2014 | Volume 8 | Article 91 | 9 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 | Mean estimated SNR N T /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.
progressively shorter bin-widths as the speed of walking increases. During standing, the peak SNR N T /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).
In contrast, the neural activity from other units in the leg area are characterized by SNR N T /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 SNR N T /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 SNR T /T toward longer rather than shorter bin-widths when the Monkey is running compared to walking (Figure 9C). The convex shape of the SNR N T /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 Figure 9. PSDs were computed from the binned neural activity with 100 ms bins and frequency resolution of 0.039 Hz.

FIGURE 10 | Power Spectrum Density (PSD) of the neural activity analyzed in
SNR N T /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 SNR N T /T and mean normalized SNR N T /T across significantly modulated units with peak above 0.5 s −1 (n = 48) depict a wide peak around 50-110 ms, in most cases. The only exception is the normalized SNR N T /T during standing, which is highest at the shortest bin-width analyzed (25 ms). Around the peak, the magnitude of the SNR N T /T are similar independent of the speed so the curves overlap each other. At longer bin-widths SNR N T /T during slow walking is higher than during fast walking or running, and is highest during standing.

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 binwidths. 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 SNR N T /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 SNR N T /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 SNR N T 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 SNR N T /T curve. Thus, the dead-time limits the benefit of reducing the bin-width, since the estimated SNR N T becomes increasingly distorted. The bin-width at which the peak occurs maximizes the trade-off between SNR N T 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 SNR N T /T curves that are also characterized by a peak. So the timescales are similarly characterized by the bin-widths at which the SNR N T /T curves peak. In contrast, when the recovery period is characterized by enhanced, rather than inhibited, probability of firing, the SNR N T /T over-estimates (rather than under-estimates) the SNR T /T. Thus, as the bin-width decreases, SNR N T /T increases without peaking or saturating. In those cases, the appropriate time-scale is not conclusive from the proposed analysis, since the higher SNR N T /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 SNR N T /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 SNR N T /T curve remains the same.
The observation that RM-SNR increases when switching to brain control indicates that the variance of the rate increasesand 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 SNR N T /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 binwidth 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 SNR N T /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 SNR N T /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 SNR N T /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 deadtime 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 SNR N T /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 SNR N T /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.
Substituting Equations (B1) and (B3) in Equation (4), the estimated SNR from a dead-time modified DSPP is: 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): For τ dλ (σ ) 1, the integral defining¯ T can be approximated by:¯ T ∼ = 1 1+τ d λ 0 t+T/2 t−T/2 λ(σ )dσ = T 1 + τ d λ 0 , so its statistics can be approximated by: Hence, the SNR T can be related to SNR T by: Using a similar approximation, the second term in Equations (B2-B4) can be approximated as: 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: Note that the second term increases with the dead-time both absolutely and as a fraction of the first term.