Abstract
The successful development of motor neuroprosthetic devices hinges on the ability to accurately and reliably decode signals from the brain. Motor neuroprostheses are widely investigated in behaving non-human primates, but technical constraints have limited progress in optimizing performance. In particular, the organization of movement-related neuronal activity across cortical layers remains poorly understood due, in part, to the widespread use of fixed-geometry multielectrode arrays. In this study, we use chronically implanted multielectrode arrays with individually movable electrodes to examine how the encoding of movement goals depends on cortical depth. In a series of recordings spanning several months, we varied the depth of each electrode in the prearcuate gyrus of frontal cortex in two monkeys as they performed memory-guided eye movements. We decode eye movement goals from local field potentials (LFPs) and multiunit spiking activity recorded across a range of depths up to 3 mm from the cortical surface. We show that both LFP and multiunit signals yield the highest decoding performance at superficial sites, within 0.5 mm of the cortical surface, while performance degrades substantially at sites deeper than 1 mm. We also analyze performance by varying bandpass filtering characteristics and simulating changes in microelectrode array channel count and density. The results indicate that the performance of LFP-based neuroprostheses strongly depends on recording configuration and that recording depth is a critical parameter limiting system performance.
Introduction
A central challenge for the treatment of a wide range of brain disorders involves developing sensors of brain activity that can help to restore lost function (Tresco and Winslow, 2011). Brain-machine interfaces (BMIs) use implanted electrode arrays to restore lost motor function by decoding motor intention (Schwartz et al., 2006; Hatsopoulos and Donoghue, 2009; Scherberger, 2009; Andersen et al., 2010; Gilja et al., 2011). Presently, the best-performing BMIs are based on neuronal spike recordings (Hochberg et al., 2006; Santhanam et al., 2006; Velliste et al., 2008; Aggarwal et al., 2009; Ganguly et al., 2011; Simeral et al., 2011), however, local field potentials (LFPs) offer promise as a robust, reliable alternative to spike signals (Pesaran et al., 2002; Andersen et al., 2004; Scherberger et al., 2005; Zhuang et al., 2010; Bansal et al., 2011). Unlike recordings of action potentials, which are generated by individual neurons, LFPs are a continuous, pooled measure of neural activity in the vicinity of the recording electrode (Xing et al., 2009) and are present in the absence of action potential recordings (Bokil et al., 2006). Therefore, LFP activity has the potential to offer greater signal reliability than spike recordings.
Current LFP-based systems have not been optimized for LFP recordings, and it is unknown whether or how the performance of these systems is limited by the design specifications of chronically implanted microelectrode arrays. For example, although it is known that the properties of LFP recordings depend on the depth of the recording electrode within cortex (Xing et al., 2009; Maier et al., 2010), the relationship between recording depth and the performance of LFP-based BMIs has not been previously investigated. The use of fixed geometry microelectrode arrays, rather than movable microelectrode arrays, discourages such investigation because optimizing the properties of recordings involves testing each recording configuration with a new implant. As a result, optimizing LFP-based neural sensors has been impractical.
Here, we use a chronically implanted microelectrode array with 32 independently movable electrodes to optimize the decoding of movement goals from neural activity (Gray et al., 2006; Gray and Goodell, 2007, 2009). We recorded spiking and LFP activity from prefrontal cortex while two monkeys performed a memory-guided saccade task. Over a span of months, we systematically advanced the depth of the recording electrodes through the prearcuate gyrus by up to 100 μm each day in both animals. For each recording session, we decoded the goal of the saccadic eye movement from neuronal activity during the memory period before the saccade. The results indicate that recording depth is a critical parameter that constrains the performance of LFP-based neural interfaces. We present strategies for improving performance by optimizing recording depth, bandpass filtering characteristics, and microelectrode array properties such as channel count and density. A preliminary report of these data was presented in (Markowitz et al., 2011).
Materials and Methods
Experimental preparation
Two adult male rhesus macaques (Macaca mulatta) participated in the study (Monkey A and Monkey S, 9.5 kg and 8.4 kg, respectively at the start of the experiments). Both animals had been used previously in other experiments studying eye movements, and identical training protocols were used for both animals (see Memory-guided saccade task). Before behavioral training, each animal was instrumented with a head restraint prosthesis to allow fixation of head position and tracking of eye position. Each monkey was behaviorally trained for several weeks in an unlit sound-attenuated room (ETS-Lindgren). Following behavioral training, we implanted a low-profile recording chamber (Gray Matter Research) in a craniotomy made over the right prearcuate cortex of each animal using image-guided stereotaxic surgical techniques (Brainsight, Rogue Research). We also installed a stainless steel ground screw over visual cortex to provide a non-local point of reference. A semichronic microelectrode array microdrive (SC32-1, Gray Matter Research) was then inserted into the recording chamber and sealed (see Microelectrode microdrive design). All surgical and animal care procedures were approved by the New York University Animal Care and Use Committee and were performed in accordance with the National Institute of Health guidelines for care and use of laboratory animals.
Memory-guided saccade task
Each monkey performed a memory-guided saccade to one of eight targets for a liquid reward. All trials began with the illumination of a central red square, on which the animal needed to fixate for a baseline period (500–800 ms). A second red square, the “cue,” was flashed for 300 ms after the baseline period to a peripheral location to indicate the target of the saccade. After a delay period (0.7–1.2 s for Monkey A, 1–1.5 s for Monkey S), the central fixation square was extinguished, providing the go signal for the animal to saccade to the peripheral cue location. Within 100–150 ms after the saccade, the cue reappeared and the animal had to fixate on the red square for an additional 300 ms. On each trial, the target was randomly presented at one location on a grid of eight possible locations spaced 10° around the central red square. Target locations were interleaved trial-by-trial in equal proportions. During each behavioral session, Monkey A performed 120–150 memory-guided saccades, and Monkey S performed 200–250 saccades.
A trial was aborted if the monkey failed to align its gaze within 2° of the center of the fixation or saccade targets. When an abort was detected, all visual stimuli were extinguished immediately, no reinforcers were delivered, and the trial was restarted after a 1200–1800 ms intertrial interval. Both monkeys rarely aborted trials (4% for Monkey A, 5% for Monkey S). Aborted trials were excluded from further analyses. Data reported here were collected after at least 3 weeks of training on the memory-guided saccade task.
Microelectrode microdrive design
All data reported here were obtained using a semichronic microelectrode array microdrive, SC32-1 (Gray Matter Research). The SC32-1 is a micromanipulator system capable of independent bidirectional control of 32 microelectrodes. The system is designed to be chronically implanted within a recording chamber system (Fig. 1a). Electrodes are spaced by 1.5 mm and each has a travel range up to 20 mm. Each actuator consists of a lead screw fixed within a housing and an eccentric brass shuttle mounted to the electrode (Fig. 1b). A complete turn of the lead screw advances the electrode 125 μm.
Recording protocol
At the time of implantation, the initial position of each recording electrode was recessed ∼1 mm within the drive. Electrodes were advanced through a Silastic membrane in the recording chamber, the dura mater and pia before entering the cortex. Action potentials were first recorded at a median depth of ∼2–3 mm beyond their initial position (2.23 mm in Monkey A; 3.04 mm in Monkey S). Electrodes were gradually advanced between daily recording sessions (mean 34 μm/d in Monkey A; 100 μm/d in Monkey S) until action potentials were no longer present, indicating passage into white matter. We obtained neural recordings by advancing the electrodes up to a median distance of 6 mm from their initial position. Neural recordings during the performance of the behavioral task were first obtained in Monkey A when clear multiunit activity was present on all channels, i.e., at the top of cortex. In Monkey S, neural recordings during behavioral task performance were initially obtained at more superficial depths above the cortical surface, before multiunit activity was present.
Data acquisition
Eye position was continuously monitored with an infrared optical eye tracking system sampling at 120 Hz (ISCAN). Eye positions were digitized at 1 kHz. Visual stimuli were presented on an LCD screen (Dell Inc.) placed 34 cm from the subjects' eyes. The visual stimuli were controlled via custom LabVIEW (National Instruments) software executed on a real-time embedded system (NI PXI-8184, National Instruments). Neural recordings were made with glass-coated tungsten electrodes (Alpha Omega) with impedance 0.7–1.5 MΩ measured at 1 kHz (Bak Electronics). Neural signals were preamplified (10× gain; Multi Channel Systems), amplified and digitized (16 bits at 30 kHz; NSpike, Harvard Instrumentation Laboratory, LSB = 1 μV after preamplification), and continuously streamed to disk during the experiment (custom C and Matlab code).
LFP waveforms were computed from the broad-band activity by low-pass filtering the raw, broad-band recording at 300 Hz and down-sampling at 1 kHz. Multiunit firing rates were obtained by high-pass filtering the raw data at 300 Hz, maintaining the original 30 kHz sampling rate, and identifying threshold crossing events, as described below. In later analyses, multiunit activity (MUA) is treated as a continuous high-frequency voltage signal to enable direct comparison with LFP data.
Data analysis
Depth registration.
To characterize decoding performance when all electrodes are at corresponding cortical depths, we aligned the depth of recordings across sessions. After the completion of all experiments and all 32 electrodes had been advanced into white matter, we registered previously measured absolute cortical depths to a common zero point across the array (i.e., the cortical surface) using an iterative optimization algorithm. First, for each channel, we estimated the variance of the LFP signal from a 60 s recording obtained at each sampled depth while the animal was sitting quietly. Once a complete profile of LFP variance versus absolute depth was obtained, we performed a variance stabilizing transformation to correct for minor differences in scale due to variable electrode impedances across the array. This transformation involved calculating the logarithm of the LFP variance at each depth and then normalizing log variance data to range from 0 to 1 on each channel. Next, we identified the alignment for each channel pair that minimized the pairwise Euclidean distance between their normalized variance versus depth profiles. To find this optimal alignment, we calculated the pairwise Euclidean distance at vertical offsets ranging from −1 mm to +1 mm. The offset with minimal pairwise Euclidean distance was labeled “optimal.” To determine the best vertical offset for a given channel, we estimated the optimal offset with respect to each of the other 31 channels, and then calculated the mean of these values. Then we shifted the variance versus depth profile for each channel by half of its optimal mean offset. This shift operation was performed in batch mode, with all electrodes shifted simultaneously once all mean offsets had been calculated. This entire procedure was repeated iteratively until no further shifts were required. Results were inspected by calculating the array-averaged MUA and LFP variance versus depth profiles and looking for inflection points characteristic of sudden changes in activity at specific depths. We labeled this inflection point zero cortical depth in both animals, and confirmed its correspondence to the top of cortex by plotting the depths of all observed spiking units across the array. In both animals, the inflection point occurred just above the first observed spikes on all channels.
Spectrographic analysis.
We estimated LFP power as a function of frequency and time using multitaper spectral estimation (Mitra and Pesaran, 1999) with 15 tapers, 20 Hz smoothing, and a 400 ms estimation window. We selected this degree of smoothing to give reliable estimates of LFP power on individual trials. When we reduced the degree of smoothing, the increased variability in power estimates reduced decoding performance.
Multiunit rate analysis.
To identify multiunit spike events, we first bandpass filtered the broad-band continuous raw recording from 300 to 5000 Hz. We then applied a 3.5 SD threshold to detect putative spike events. All waveform peaks that exceeded this threshold were then labeled as multiunit spike events, and the corresponding timestamp was recorded for later analysis. We estimated multiunit rate by counting the number of multiunit events during a 400 ms decoding window immediately after target offset and multiplying by 2.5. This provided a rate estimate in events/s. The duration of the decoding window was selected to match the analysis of LFP activity and enable a direct comparison of both signals. Similar decoding results were obtained when we varied the threshold for detection up to 4.5 SDs.
LFP decoding.
We decoded LFP array recordings from single trials to predict the saccade movement goal from one of eight potential locations. The decoding algorithm contained three major stages: preprocessing, feature extraction and linear decoding (Fig. 2). During preprocessing (Fig. 2a), LFP traces from each trial are transformed from the time domain to the frequency domain and reshaped into an Nch × Nf × 1-element row vector, Xtrial. This sample vector is then projected onto a Template matrix, M, during the Feature Extraction stage (Fig. 2b). The Template contains the first N modes (typically, 30) of a previously analyzed training dataset, and serves the function of reducing the dimensionality of Xtrial to a set of features that capture most of the variance in the original high-dimensional space.
The scatter plot shown in Figure 2b illustrates how artificially generated sample data from eight target classes (represented by different colors) might be projected onto a two-dimensional subspace. This projection reveals tightly clustered samples from a few of the classes, which might be useful for decoding purposes. The poorly clustered samples may not be distinguishable in this subspace, but further projection of the training data onto additional modes may reveal clear boundaries in a higher dimensional space.
Features were calculated from memory delay activity during a 400 ms window beginning immediately after the target cue was extinguished. To limit the influence of the evoked response due to target offset, we subtracted the mean waveform across trials from the activity before subsequent spectral analysis. Only data during the memory period, after target offset and before the go cue, was included in the decoding procedure. The log power spectrum from 0–300 Hz was estimated on each of the 32 electrodes using multitaper methods, as described above. The LFP time series was padded to give 327 frequency bins yielding a total of 32 × 327 = 10,464 features. We then used singular value decomposition (SVD) to linearly project the signal onto a 30 dimensional subspace that was used for classification analysis (see below). In subsequent analyses, we explored changes in decoding performance using a variable maximum passband from 10 Hz to 10 kHz and a range of subspace dimensions from 5 to 80 modes.
The 30 dimensional subspace of SVD modes was calculated from a training dataset containing repeated interleaved presentations of all eight saccade goals. We performed a spectral analysis of the delay activity during all training trials, yielding a Ntrials × 10,464 dimensional rectangular matrix, Xtrain. We then calculated the inner product between this training matrix and itself, yielding a square training matrix, XXT, with dimension Ntrials × Ntrials. Finally, we performed a singular value decomposition of the square training matrix: The spectral modes of the training data were calculated by projecting the transpose of the rectangular training matrix onto the first 30 right singular vectors of the SVD: Linear Decoding occurs at the final stage of the algorithm (Fig. 2c). Decoding of the target location on each trial proceeds by calculating a 10,464 element spectral array, Xtrial, from activity during the memory epoch (300–700 ms after cue onset), projecting this array onto the 30 spectral modes described above, and performing a linear discriminant analysis (Duda et al., 2000) to classify which target is most likely given the test data: where wi is a projection matrix that is used to optimally discriminate whether or not Xtrial was drawn from the “Target i” class, and the function f(x) is a nonlinear mapping from the discrimination subspace onto a probability value for each predicted location, which we treat as the Probability of Decoding Targeti from a trial sample. After estimating this probability for all eight targets, an argmax operation is applied to identify the most likely decoding classification. The decoded target direction is then used to predict where the monkey is planning to move his eyes. We used the “classify” command in Matlab to construct a simple linear decoder from the training data and a corresponding array of saccade target labels. Classifier performance estimates were bootstrapped using leave-one-out cross-validation. Model performance during each experimental session was summarized by the mean correct performance averaged across all movement goals, and by a confusion matrix quantifying the probability of predicted target directions, conditioned on all observations within each target class.
LFP decoding by spectral band.
To decode movement plans for specific frequency bands, we calculated the mean LFP power in the spectral range of interest on each channel, yielding 32 features on each trial. Then we used SVD to identify the modes of this reduced-dimensionality dataset before applying the previously described decoding algorithm. Typically, maximum performance was achieved using five modes. It is important to note that these modes reflect spatial patterns of activity across the 32-channel array in a limited spectral band, rather than high-dimensional structure in a 10,646-dimensional channel-frequency feature space.
Multiunit rate decoding.
To decode movement plans from multiunit firing rate estimates, we used data samples with 32 features, representing the multiunit firing rate observed on each electrode during a given memory epoch. This reduced-dimensionality information was then used in place of the 10,464-dimensional LFP data in the linear decoding procedure described above.
Decoding at registered depths.
To study decoding performance at corresponding cortical depths across the array, we developed a procedure for constructing “virtual” sessions from discontinuously recorded data. After choosing a specific registered depth for study, we identified the session during which each electrode was closest to this location and selected the corresponding neural data from that channel and recording day. Typically, neural data were drawn from 5 to 10 unique sessions, and all selected channel data were <200 μm from the target depth. Finally, we grouped voltage traces from all 32 channels to create virtual trials, such that all 32 traces assigned to a given trial were associated with the same cue location in their original recording sessions. Throughout this study, we use the term “registered cortical depth” when describing virtual session data, and “mean electrode depth” to describe the mean absolute depth of electrodes in simultaneously recorded data. It is important to note that both of these terms refer to the depth in cortical tissue and may not reliably correspond to depth within the cortical sheet. Although the microdrive was implanted approximately normal to the gyral surface in both animals, some electrodes may have penetrated sulcal banks and remained in the same cortical layer over a span of several millimeters.
N-channel performance estimation.
We studied the influence of channel count, Nchannels, on decoding performance by randomly selecting subsets of channels from the same experimental session when the analysis called for Nchannels < 32. When the analysis called for Nchannels > 32 we pooled channel data from consecutive experimental sessions. Decoding performance reported for Nchannels < 32 data are averages over classifiers constructed from 20 randomly selected subsets of channels. Reported data are the maximum performance observed by building decoders using from 5 to 80 modes of the training dataset. In general, it is necessary to increase the number of degrees of freedom to improve performance when Nchannels > 32.
Pairwise correlation analysis.
We investigated the relationship between electrode separation and decoding performance using pairwise LFP correlation as a linking variable. This approach is based on the intuition that pairwise correlation should approach zero as electrodes approach infinite separation and one as they converge to a single point in space. By extension, we expect decoding performance to drop as electrodes converge and the information content drops to that of a single location. We incorporated these assumptions into the two-step procedure described in Equations 4–7 below. In the first step of our analysis, we linearly combined 32 channels of LFP data from a single session to yield the mean LFP correlation level that would be expected at a specific mean electrode separation. After mean-subtracting and normalizing all 32 LFP traces to have the same SD, we chose a reference channel from the array as the physical point in space toward which all other electrodes would be virtually shifted. Then we applied the following series of “blending” transformations to the remaining 31 “test” channels to yield the increased mean pairwise correlation that would be expected as all electrodes converged toward the reference point: where ρ is the desired mean pairwise correlation and ζi(t) is a Gaussian-distributed stochastic process with property <ζi(t) · ζi(t)T〉 = 0 that is generated randomly for each test channel. As the value of ρ increases, Equation 4 destroys channel-specific information while Equation 5 reshapes the test signal to look more like the reference signal. Equation 4 incorporates the biophysical assumption that electrical potentials decay in amplitude and are increasingly corrupted by noise with increasing distance from a point source. Put differently, as a test electrode is moved closer to the reference point and further away from its original location, the component of the LFP that derives from the original recorded waveform should look progressively weaker and noisier. Equation 5 incorporates our assumption that the test signal should look more like the reference waveform as it approaches this location.
For each choice of reference channel, this algorithm was applied to all trial data for the other 31 test channels. Then we trained a linear decoder on the blended data to quantify decoding performance at the desired correlation level. This procedure was repeated iteratively, using each channel as the reference, to limit bias from any one cortical site in our final results. Each decoding performance value reported in the text is an average over all 32 runs of this analysis. Each pairwise correlation value reported in the text was calculated directly from the blended data.
After calculating mean decoding performance versus mean pairwise correlation, the second step of our procedure used a simple model to map correlation values onto electrode distances. This model was motivated by experimental measurements of this relationship by (Destexhe et al., 1999) and more recently by (Eggermont et al., 2011), which suggest the following relationship: The parameter of this model is determined by the mean pairwise correlation observed in our original untransformed data and the fixed 1.5 mm distance between adjacent electrode pairs. We used this model to map pairwise correlation values in our blended data onto estimated electrode distances. This enabled us to plot mean decoding performance against estimated distance.
Results
To optimize the performance of a neural decoder for neuroprosthetic control, it is important to understand the functional organization of movement plans as a function of cortical depth. We investigated this problem by training two adult male rhesus macaques (Macacca mulatta) to perform a memory-guided saccade task (see Materials and Methods). Both animals were then instrumented with a circular recording chamber system and implanted with a chronic 32-channel microdrive over prearcuate cortex (Fig. 1). The span of the multielectrode array covered the frontal eye fields (FEF) and dorsolateral prefrontal cortex (dlPFC) in both animals.
Semichronic microdrive recording stability
We recorded 42 behavioral sessions over a span of 59 post-implantation days in Monkey A, and 56 sessions over 68 d in Monkey S. The microdrive provided stable recordings of MUA and LFP signals over the entire post-implantation interval in both animals (Fig. 3). Changes in the signal variance across consecutive recording days typically reflected the advancement of electrodes into more active zones of cortex. We did not observe a sustained influence of post-implantation time on the quality of spike recordings (Fig. 3aii) or MUA variance in either animal (Fig. 3b). Normalized variance was used to control for differences in signal variance across electrodes, for example due to differences in electrode impedance. For Monkey A, average normalized MUA variance did not change from the beginning to the end of the experiment (0.49 ± 0.01 SEM during the first and last 20 d), while average normalized LFP variance showed a moderate decrease of 0.03 (0.48 ± 0.01 SEM during the first 20 d vs 0.45 ± 0.01 SEM during the last 20 d). For Monkey S, average normalized MUA variance increased by 0.01 from the beginning to the end of the experiment (0.46 ± 0.01 SEM vs 0.47 ± 0.01 SEM), while average normalized LFP variance increased by 0.11 (0.33 ± 0.01 SEM vs 0.44 ± 0.03 SEM). The sustained increase in LFP variance after day 40 in Monkey S corresponds to mean electrode depths >2 mm, likely outside the gray matter of the superficial prearcuate gyrus.
Depth dependence of LFP decoding performance
We mapped the functional organization of movement plans as a function of cortical depth by decoding 32-channels of simultaneously recorded neural activity from each behavioral session. We used a linear, multivariate spectral decoding algorithm to decode movement direction (1 of 8 possible directions) from single trials (Fig. 2; see Materials and Methods). We first present the performance of our decoding algorithm by generating confusion matrices for superficially recorded sessions in both animals (Fig. 4a; Monkey A: 88 trials, median depth = 0.126 mm; Monkey S: 247 trials, median depth = 0.408 mm). These matrices plot the probability of predicted target directions, conditioned on all observations within each target class. Cells along the diagonal indicate the probability of correct predictions, while off-diagonal cells indicate errors. In both monkeys, decoding errors typically predicted directions that were similar to the actual direction. Overall decoding performance showed a strong bias toward contralateral (leftward) directions in both animals (Fig. 4b), consistent with previously observed workspace effects in macaque frontal cortex (Funahashi et al., 1989). Mean decoding performance across all targets is 83% for both animals.
The solid black traces in Figure 4c show mean decoding performance across all eight targets for each recording session in both monkeys. These data are reported as a function of mean absolute electrode depth across 32 channels of simultaneously recorded channel data, where zero represents the site at which spiking activity was first observed on each channel. Therefore, negative depths correspond to locations above the cortical surface. The extent of superficial depths studied differed between animals due to differences in the recording protocol (see Materials and Methods). Peak decoding performance is 83% in both animals, and the corresponding “best” sessions are shown in the preceding panels (Fig. 4b). Since our decoding algorithm achieves these results by operating on the power magnitude at each LFP frequency, and otherwise discards timing information contained in the phase, we reasoned that performance would be improved or at least maintained by eliminating trial-specific correlations in LFP power across the array. To study this question, for each electrode, we shuffled LFP traces across trials within each target class. This procedure maintains the mean power across the array for each target while destroying correlations in power between electrode pairs across trials. The dotted traces in Figure 4c show mean decoding performance using shuffled data. These results are generally comparable to performance achieved using unshuffled data, shown by the solid traces in Figure 4c. Decoding performance using shuffled data is significantly less than performance using unshuffled data at 5% of sampled depths in Monkey A (n = 40 depths ranging from −0.2 to +1.5 mm, p < 0.05) and 33% of sampled depths in Monkey S (n = 45 depths ranging from −1 to +3 mm, p < 0.05). The average decrease in performance across all depths by shuffling is −0.05 ± 0.01 SEM in Monkey A and −0.11 ± 0.02 SEM in Monkey S. These results demonstrate that LFP power estimates can be drawn from different trials with only a 5–10% expected decrease in performance. The results also suggest that correlations between LFP power on different channels in these data do not impact performance, which we will examine in more detail below.
LFP decoding varies with recording depth more than multiunit decoding
The data in Figure 4c hint at a trend toward decreasing performance with lower depth. However, since different electrodes in the array typically were advanced by different amounts each day, the dependence of decoding performance on mean electrode depth across the array is potentially misleading. This is because not all electrodes were registered to the same offset from zero cortical depth during each recording session. To study the representation of movement goals at corresponding depths across the array, we resampled data from each animal to create virtual sessions. Each virtual session was constructed by grouping 32 LFP traces that were recorded at the same cortical depth and for the same target at all locations in the array, but which may have been drawn from different recording sessions. This approach can destroy correlations in activity by pooling discontinuously recorded channel data, and therefore is comparable to the previously reported shuffled analysis.
We studied the mean decoding performance at registered cortical depths in both monkeys (Fig. 5a) for the following target categories: all eight targets (black lines), three targets in the contralateral hemifield (red lines), and three targets in the ipsilateral hemifield plus two on the vertical axis (blue lines). We present decoding performance for both hemifields to quantify the strength of workspace effects in this region of the brain. In general, this is important for understanding the factors limiting overall decoding performance, and can help guide the design of more effective recording and decoding strategies. In Monkey A, peak performance across all eight targets is 80% at −0.10 mm and minimum performance is 28% at 1.25 mm (Fig. 5ai). The data in this panel show a substantial decrease in performance at deeper cortical sites. In Monkey S, we explored a broader range of depths, from −1 mm above the cortical surface to +3 mm into the white matter (Fig. 5aii), and recovered similar results: peak performance across all eight targets is 71% at 0.50 mm and minimum performance is 15% at −0.70 mm above the zero point, and 20% at 2.90 mm below the zero point. The data from Monkey S show relatively poor LFP decoding performance at surface cortical locations, which increases as electrodes enter superficial cortex and then decreases at deeper cortical sites. In both panels, decoding performance is typically much better for contralateral targets than ipsilateral targets. Averaged across all registered depths, contralateral LFP decoding performance is 0.20 ± 0.12 SEM (n = 60 depths) higher than ipsilateral performance in Monkey A and 0.09 ± 0.10 SEM (n = 60 depths) higher in Monkey S.
To establish a quantitative link between LFP and MUA decoding performance at corresponding depths, we repeated the above analysis for MUA by adapting our decoding algorithm to operate on a feature vector of 32 multiunit firing rate estimates on each trial (Fig. 5b). In Monkey A (Fig. 5bi), peak MUA decoding performance across all eight targets is 80% at −0.10 mm and minimum performance is 28% at 1.25 mm. The data in this panel show a decrease in performance at deeper cortical sites, similar to the trend in Figure 5ai. In Monkey S (Fig. 5bii), peak MUA decoding performance across all eight targets is 71% at 0.50 mm and minimum performance is 15% at −0.70 mm above the zero point, and 20% at 2.90 mm below the zero point. These data show the same general trend observed in Figure 5aii, including the relatively poor MUA decoding performance at surface cortical locations, which increases as electrodes enter superficial cortex and then decreases at deeper cortical sites. Averaged across all registered depths, contralateral MUA decoding performance is 0.19 ± 0.06 SEM (n = 60 depths) higher than ipsilateral performance in Monkey A and 0.13 ± 0.13 SEM (n = 60 depths) higher in Monkey S.
The above analysis demonstrates that LFPs and MUA are most informative about eye movement goals and yield comparable decoding performance at superficial cortical depths. To understand the relative decoding performance of LFP and MUA activity at superficial and deep cortical sites, we performed a spectral decoding analysis that systematically varies the maximum passband frequency, Fmax (Fig. 6ai,bi). In contrast to the preceding multiunit rate analysis, this spectral analysis treats the MUA signal as a continuous voltage signal, rather than a sequence of discretely thresholded events. At superficial sites in both animals, decoding performance reaches its maximum value when Fmax is increased to 300 Hz (the upper bound used here to filter LFP signals from raw data). Interestingly, inclusion of multiunit power at frequencies up to 1 kHz does not yield any performance improvement. Similar results for multiunit signals from 1 kHz to as high as 10 kHz were also observed (data not shown). This indicates that the information conveyed by MUA through high-frequency oscillations is consistent with LFP activity when decoding eye movement goals from activity at superficial cortical depths.
The performance similarity between LFP and MUA is greatest at superficial depths and is not present at deep sites. At deep sites (∼1 mm below the superficial sites in both animals), decoding performance continues to improve as Fmax increases >300 Hz. To enable a more direct comparison of these superficial and deep sites, we rescaled the corresponding 8-target mean decoding performance profiles to range from 0 to 1 (Fig. 6aii,bii). The vertical distance between these rescaled traces indicates that LFPs contain more information about movement goals at superficial sites than at deep sites. However, to maximize decoding performance at deep cortical sites, a larger passband that recovers multiunit spiking is required. Since LFP decoding performance decreases more substantially than multiunit decoding performance at deeper sites, it is unlikely that LFP activity above ∼100 Hz strongly reflects multiunit activity.
Variations in LFP decoding with electrode spacing and electrode number
The spacing of electrodes in an array is likely to influence the correlation structure of recorded neuronal activity, which may in turn constrain decoding performance. To understand how LFP decoding performance is influenced by correlated activity across the array, we examined how 8-target decoding performance in the best performing session varies as a function of the spacing between electrodes in the recording array (Fig. 7; see Materials and Methods). We began by calculating the mean pairwise LFP correlation across the array at 1.5 mm spacing, and then reshaped LFP waveforms to mimic the correlation structure that would be expected as electrodes are drawn progressively closer together. Using simultaneously recorded data from the optimal recording session in each animal, we artificially increased these correlations in sequential rounds and calculated decoding performance in each case (Fig. 7a). As expected, when LFP correlations approached 1, LFP decoding performance approached the performance obtained with a single electrode.
We established a quantitative connection between decoding performance and electrode spacing using a simple model that maps pairwise LFP correlation onto estimated recording electrode separation (see Materials and Methods). By plotting mean decoding performance versus estimated electrode spacing in mm, we were able to identify the spacing at which performance decays to 95% of the value observed at 1.5 mm (Fig. 7b): 0.35 mm in Monkey A and 1 mm in Monkey S. These results suggest that decoding performance would not be strongly affected by electrode spacing slightly <1.5 mm. However, when 32-channel electrode spacing goes <0.5–1 mm, 8-target LFP decoding performance falls significantly.
We reasoned that increasing the number of channels might help to increase performance in the presence of strong pairwise correlation structure. To study this problem, we varied the number of channels in the analysis from 4 to 128 by selecting channel data from the best-performing session (when Nchannels ≤ 32), or by pooling channel data across this and preceding/following sessions when Nchannels > 32). We observed that increasing the number of channels significantly increases performance in both animals (Fig. 7c). In both animals, performance saturates at a maximum value as the number of channels is increased from 32 to 128. Although the performance improvement above 32 channels is negligible for Monkey A, Monkey S continues to show improved decoding performance for ipsilateral targets as the number of channels increases to 128.
Spectral band contributions to decoding performance
Our spectral decoder was designed to extract the most relevant information for decoding from a broadly defined frequency spectrum, without making a priori assumptions about which frequencies are most informative. However, to provide a more direct connection between our analysis and previous work, we also studied the information content of commonly studied spectral bands, such as the α (6–12 Hz), β (12–30 Hz) and gamma (30–100 Hz) bands. For each spectral band, we repeated our decoding analysis using the mean power across the corresponding frequency range on each channel (Fig. 8). Therefore, training data were reduced to a 32 × 1 dimensional array for each trial. Following this averaging procedure, we applied all steps of the algorithm described previously, and calculated the maximum decoding performance using from 5 to 30 SVD modes to define the feature subspace. Typically, optimal performance was observed using only 5 modes. This analysis reveals that spectral bands <30 Hz yield <55% mean performance across all eight targets at nearly all depths in Monkey A, and <25% performance in Monkey S (Fig. 8ai,bi). By contrast, the gamma frequency band yields consistently better performance, most notably in the 50–100 Hz band Fig. 8aii,bii), where mean performance reaches 74% in Monkey A and 50% in Monkey S. To directly connect this analysis with the data in Figure 6, we also studied decoding performance versus depth in higher frequency bands, including 100–300 Hz, 0.3–1.0 kHz, and 1.0–3.0 kHz (Fig. 8aiii,biii). As observed previously, the mean performance at depths ∼0 mm are similar in the high LFP (100–300 Hz) and MUA (0.3–1.0 kHz) bands, but substantially diverge at deeper sites. Notably, the mean performance using each spectral band at each registered depth was not consistently higher than the performance observed using our earlier analysis (e.g., peak performance in 100–300 Hz band: 74% in Monkey A, 73% in Monkey S). This suggests that no single frequency band provides all of the useful information for decoding, but rather, the information is distributed across a range of frequencies.
Discussion
Here we use chronically implanted movable microelectrode arrays to decode saccade movement goals from LFP activity recorded at different depths within the prearcuate cortex. Decoding performance is best at superficial depths, proximal to the first observed spiking activity after implantation, and LFPs and multiunit firing rates yield comparable decoding performance at these superficial sites. At deeper cortical sites, performance drops for both LFP and multiunit activity, and the reduction in performance is significantly greater for LFP activity than for multiunit firing rate. Instead of defining specific frequency bands for decoding, we used a spectral decoding algorithm that generates features using a weighted average across all frequencies. By varying filter settings from 50 Hz up to 1000 Hz, we find that at superficial sites LFP activity up to ∼300 Hz can be decoded to capture almost all of the information that is available from both LFP and multiunit activity. Importantly, however, LFP decoding performance at deeper sites is significantly worse than decoding the full-bandwidth signal. Finally, we use a model-based approach to derive a minimum electrode spacing constraint for efficient LFP acquisition. It should be noted that we do not observe a globally optimal “peak” performance within any of these parameter regimes as, for example, peak decoding performance in depth occurs over a range of 100–200 μm in both animals.
Previous work has shown that movement selectivity in LFP activity may be useful for neuroprosthetic control applications. However, the influence of cortical recording depth on system performance was previously unknown. Earlier work decoded LFP activity using acutely inserted electrodes (Pesaran et al., 2002; Mehring et al., 2003; Rickert et al., 2005; Scherberger et al., 2005; Bokil et al., 2006). More recent work has decoded LFP activity recorded on chronically implanted, fixed geometry electrode arrays (Mollazadeh et al., 2008; Zhuang et al., 2010; Bansal et al., 2011). Due to the sampling limitations of acute recordings and physical constraints imposed by fixed geometry arrays, these studies were unable to systematically quantify the relationship between decoding performance and recording depth, as we do here using movable implanted electrodes.
Comparison between LFP and multiunit activity
The multiunit decoding performance we report compares favorably with other work. We show that both LFP- and MUA-based decoders are capable of predicting eye movement goals on a single trial basis with comparable accuracy (1-of-8 performance: LFP power = 70%, 80%; Multiunit rate = 71%, 80%) from 32 channels of data. Two studies have also presented results decoding 1-of-8 movement goals during instructed-delay periods, permitting comparisons. Scherberger et al. (2005) found that reach direction could be decoded correctly on 97% of trials when pooling recordings from 89 selected neurons and 81% when pooling recordings from 125 LFP sites. More recent work using 96-electrode Utah arrays in dorsal premotor cortex found that 1-of-8 movement goals were correctly decoded from spiking activity on up to 69% or 77% of trials, depending on the duration of the trial, and hence, the decoding analysis window (Santhanam et al., 2006). LFP decoding performance was not reported in this study. Although several factors differ between these experiments and our own, the multiunit and LFP decoding performance we obtain using 32 electrodes is broadly consistent with these reports.
A particularly interesting aspect of our results is the observed change in decoding performance using signals with progressively increasing bandwidths, from 50 to 1000 Hz. Varying the passband of the low-pass filter permits us to study the relative information content of LFPs and MUA as a function of recording depth. The analysis presented in Figure 6 shows that LFP activity at up to 300 Hz is sufficient to achieve near-maximum performance at superficial locations, and demonstrates that analyzing both LFPs and MUA is necessary to obtain the highest performance when recording from deeper cortical sites. The influence of recording depth is especially interesting when considering the source of LFP signals in the several hundred Hz range. There is considerable debate on the precise contributions of neural activity to extracellular potentials in the 100–500 Hz frequency range. Lower frequency LFP signals, below ∼100–200 Hz, are generally believed to reflect postsynaptic transmembrane potentials (Klee et al., 1965; Mitzdorf, 1985; Buzsáki, 2006; Okun et al., 2010). The tuning of gamma-band activity in cat V1 measured by intracellular membrane potential is closely correlated with that of spiking activity (Azouz and Gray, 2003), which itself is associated with rapid fluctuations in subthreshold membrane potential (Azouz and Gray, 2008). High-frequency LFP fluctuations also likely reflect the leakage of action potentials and spike after potentials into the LFP signal band (Ray et al., 2008; Ray and Maunsell, 2011), the degree to which depends on the amplitude of the action potentials present in the extracellular recording (Zanos et al., 2011). Since we find that the relative information content of high-frequency LFPs signals and multiunit spiking varies significantly with recording depth, however, our data suggest that there is not an obligatory relationship between nearby spiking activity and high frequency ∼100–500 Hz LFP activity. One possible explanation is that high-frequency, synaptic activity within the dendrites of pyramidal cells underlies the laminar organization of LFP activity and is tightly linked to both spiking activity and high-frequency LFP activity. The extent to which intracellular potentials reflect dendritic fluctuations remains to be tested and future work will need to examine the contribution of synaptic potentials and spiking events to LFP signals as a function of recording depth.
Comparison between LFP activity and electrocorticograms
There is increasing interest in driving a neural interface using electrocorticogram (ECoG) activity recorded at the cortical surface and epidurally (Felton et al., 2007; Schalk et al., 2007; Chao et al., 2010; Moran, 2010; Slutzky et al., 2011), but thus-far, comparisons between these signal modalities have been limited in scope. In particular, almost all studies have not directly compared the information contained in ECoG and LFP activity using the same tasks, brain regions and subjects (but see Slutzky et al., 2011). We did not record traditional ECoG signals by placing grid electrodes on the cortical surface. However, in one animal (Monkey S), we obtained LFP activity as we lowered microelectrodes through the dura, pia and into the brain. Using identical procedures across days in this animal, we found that average decoding performance at recording sites superficial to the first observed spikes was ∼30%, while performance at sites only 300–500 μm deeper in the superficial cortex was as high as 80%. These results demonstrate that penetrating the pia and inserting electrodes as little as 300–500 μm into the cortex yields substantially better neuroprosthetic control performance than placing electrodes at the cortical surface, even when microelectrodes are used. While we do not know the precise depth at which electrodes penetrated the pia and entered cortex, obtaining LFP activity in the superficial cellular regions of cortex may significantly improve decoding performance over cortical surface recordings.
Workspace effects and implications for human neuroprosthetics
Spiking activity exhibits a strong preference for contralateral visual space in macaque FEF/dlPFC (Funahashi et al., 1989), and posterior parietal cortex (Barash et al., 1991; Scherberger et al., 2005). Consistent with this, we find that LFP decoding performance in frontal cortex is nearly 100% for contralateral targets, while performance for vertical and ipsilateral targets is substantially worse (Figs. 5⇑–7). Bilateral implants may be necessary in these regions to achieve high LFP decoding performance across the entire workspace. By contrast, hemispheric workspace biases are much weaker in spiking and LFP activity in motor cortices (Georgopoulos et al., 1988; Scott et al., 2001; Cisek et al., 2003; Rickert et al., 2005), and the representation of space for eye movements is also less lateralized in humans than macaques (Kagan et al., 2010). Therefore, LFP decoding performance using a unilateral implant in human motor cortex may exceed the ipsilateral workspace performance we report here.
Caveats and limitations
Several caveats and limitations should be kept in mind before extrapolating our results to neuroprosthetic control. First, our analysis was performed off-line using activity before executed movements. Since the upcoming movement may influence performance, future work should study performance on-line without movements. Second, we studied an oculomotor/cognitive area. The depth dependence of activity in areas involved in skeletomotor movements needs to be examined. Third, we examined the encoding of goals in a two-dimensional saccadic workspace. While saccades represent a model that captures many features and challenges of this problem, studies that involve more complex reach-and-grasp movements are needed. Finally, the overall performance of neuroprosthetic control needs to be evaluated in a broader context including reliability and robustness. The reliability and robustness of spiking signals have been the focus of many reports (Williams et al., 1999; Szarowski et al., 2003; Biran et al., 2005; Polikov et al., 2005; Suner et al., 2005; Dickey et al., 2009; Simeral et al., 2011), but LFP and ECoG signals remain relatively poorly studied (but see Chao et al., 2010). Despite these limitations, this study builds substantially on previous work by demonstrating the strong influence of depth on decoding performance, and exploring strategies for optimizing performance using variable passband, channel count and density.
Footnotes
This work was supported in part by NSF Career Award BCS-0955701, NIMH Grant R01-MH087882, NIH Training Grant T32-MH19624 (D.A.M.), a Swartz Fellowship in Theoretical Neurobiology (D.A.M.), a Career Award in the Biomedical Sciences from the Burroughs Wellcome Fund (B.P.), a Watson Program Investigator Award from NYSTAR (B.P.), a McKnight Scholar Award (B.P.), and a Sloan Research Fellowship (B.P.). We thank Gerardo Moreno for surgical assistance and Roch Comeau, Stephen Frey, and Brian Hynes for custom modifications to the BrainSight system.
- Correspondence should be addressed to Dr. Bijan Pesaran, Center for Neural Science, New York University, 4 Washington Place, Room 809, New York, NY 10003. bijan{at}nyu.edu