Abstract
Deep brain stimulation (DBS) has expanded as an effective treatment for motor disorders, providing a valuable opportunity for intraoperative recording of the spiking activity of subcortical neurons. The properties of these neurons and their potential utility in neuroprosthetic applications are not completely understood. During DBS surgeries in 25 human patients with either essential tremor or Parkinson's disease, we acutely recorded the single-unit activity of 274 ventral intermediate/ventral oralis posterior motor thalamus (Vim/Vop) neurons and 123 subthalamic nucleus (STN) neurons. These subcortical neuronal ensembles (up to 23 neurons sampled simultaneously) were recorded while the patients performed a target-tracking motor task using a cursor controlled by a haptic glove. We observed that modulations in firing rate of a substantial number of neurons in both Vim/Vop and STN represented target onset, movement onset/direction, and hand tremor. Neurons in both areas exhibited rhythmic oscillations and pairwise synchrony. Notably, all tremor-associated neurons exhibited synchrony within the ensemble. The data further indicate that oscillatory (likely pathological) neurons and behaviorally tuned neurons are not distinct but rather form overlapping sets. Whereas previous studies have reported a linear relationship between power spectra of neuronal oscillations and hand tremor, we report a nonlinear relationship suggestive of complex encoding schemes. Even in the presence of this pathological activity, linear models were able to extract motor parameters from ensemble discharges. Based on these findings, we propose that chronic multielectrode recordings from Vim/Vop and STN could prove useful for further studying, monitoring, and even treating motor disorders.
Introduction
Neurosurgical implantation of deep brain stimulation (DBS) electrodes is an efficacious treatment for both Parkinson's disease (PD) and essential tremor (ET) (Ondo et al., 1998; Koller et al., 2001; Kumar et al., 2003; Rodriguez-Oroz et al., 2005; Deuschl et al., 2006). Common DBS targets include the subthalamic nucleus (STN; for PD patients) and the ventral intermediate nucleus of thalamus (Vim; for ET patients). Both structures are involved in motor control (Parent and Hazrati, 1995; Guillery and Sherman, 2002). The dorsolateral STN receives afferents from motor cortex, premotor cortex, and supplementary motor areas (Parent and Hazrati, 1995; Hamani et al., 2004). Vim projects to these areas as well as receiving afferents from the ipsilateral cerebellum.
Individual human Vim/STN neurons are active during voluntary and passive movement, somatosensation, and motor imagery (Lenz et al., 1990, 1994, 2002; Raeva et al., 1999; Magariños-Ascone et al., 2000; Magnin et al., 2000; Rodriguez-Oroz et al., 2001; Abosch et al., 2002; Benazzouz et al., 2002; Theodosopoulos et al., 2003; Williams et al., 2005). Many reports have examined single-unit activity in these regions with respect to tremor (Lenz et al., 1988, 1994, 2002; Zirh et al., 1998; Magariños-Ascone et al., 2000; Magnin et al., 2000; Rodriguez-Oroz et al., 2001; Brodkey et al., 2004; Hua and Lenz, 2005; Amtage et al., 2008) and pathological synchronous oscillations (Levy et al., 2000, 2002; Amirnovin et al., 2004). However, the number of simultaneously recorded cells in these studies was low (≤2), providing limited information about the pathological activity of larger neuronal ensembles.
Ensembles of simultaneously recorded neurons have been used to enable brain-machine interfaces (BMIs) for neuroprosthetic control in animal models (Chapin et al., 1999; Wessberg et al., 2000; Nicolelis, 2001; Serruya et al., 2002; Taylor et al., 2002; Hochberg et al., 2006; Fetz, 2007; Patil and Turner, 2008; Nicolelis and Lebedev, 2009; Lebedev et al., 2011; O'Doherty et al., 2011). However, recordings in humans have rarely been extracted from neuronal ensembles (Kennedy and Bakay, 1998; Kennedy et al., 2000, 2004; Patil et al., 2004; Quiroga et al., 2005). Human subcortical ensembles have only been used in a sole study (Patil et al., 2004), in which our laboratory demonstrated the feasibility of a subcortical motor BMI using motor control signals from neuronal ensembles recorded in motor thalamus [Vim/ventral oralis posterior motor thalamus (Vop)] or STN during DBS surgery to decode modulations of hand force during a one-dimensional target-tracking task.
We now examine neuronal population firing patterns during a voluntary motor task in a new sample of 25 human patients. While patients performed contralateral hand movements to acquire visual targets, up to 23 subcortical neurons were recorded simultaneously and acutely in Vim/Vop and STN to elucidate the relationship between neuronal modulations, rhythmic oscillations, and neuronal synchrony. A linear decoder model was applied to reconstruct cursor position from spiking activity. Neurons were classified by oscillatory firing patterns, tremor association, synchrony, and tuning to target and movement parameters.
Materials and Methods
Intraoperative recordings were conducted in 25 human patients undergoing placement of therapeutic DBS implants in either Vim or STN. All studies were approved by the Duke University Institutional Review Board and human ethics committees, and all participating patients understood and signed all required consent forms.
Patient characteristics and operative plan.
All patients selected for this study underwent either Vim (N = 14, 10 male, 4 female) or STN (N = 11, 10 male, 1 female) DBS electrode implantation surgery. Patients whose symptom presentation was dominated primarily by medication-resistant tremor (either essential tremor or severe parkinsonian tremor) were candidates for implantation in Vim, while patients with severe PD (typically akinetic/rigid variant with on/off fluctuations and dyskinesia) were candidates for implantation in STN. Both groups of patients were off their medications before and during surgery.
Patients first underwent Leksell frame placement, followed by a MRI scan to localize the implantation target. For Vim patients, the target was typically estimated according to anterior–posterior commissure (AC-PC) criteria, located ∼5–6 mm in front of the PC, on the AC-PC line with a lateral measure depending on the width of the third ventricle (typically 12–15 mm). Figure 1A shows a typical electrode trajectory. The first penetration for approaching Vim was the localization pass from a frontal burr hole, with the upper 5 mm of the recording track near the border of Vop and Vim, and the lower 5 mm in Vim and close to ventralis caudalis at the most posterior extent. Typically, the upper 5 mm was the best for multineuron recordings, reflecting more of the anterior motor thalamus (Vop rather than Vim), as shown in Figure 1A. For STN patients, the target was calculated by indirect methods, based on the AC-PC and 1 mm axial cuts of spoiled gradient recalled acquisition in the steady-state imaging. The initial target was located at 11–12 mm from the midline, 2–3 mm posterior to the midpoint of the AC-PC line, and 4 mm below the AC-PC line (Deuschl et al., 2006). Single-unit recordings were first performed to define the borders of the STN, according to standard electrophysiological criteria, with the goal of attaining at least 5.5–6 mm of STN. Typically, 2–3 passes were performed. Localization was performed using single-channel tungsten microelectrodes.
A, Sagittal diagram of the human brain at 15 mm off the midline. The vertical line represents the midpoint of the AC-PC line, which is arbitrarily set here as 24 mm in total length. The distance from the midpoint to the PC is 12 mm. The angled line represents a typical electrode trajectory, with the box emphasizing the best location (in terms of cell density) for the placement of the electrode array. This area overlaps Vim/Vop in the motor thalamus. Note that the DBS electrode would typically be placed deeper, at the border of Vim and ventralis caudatus (Vc, the sensory nucleus). The exact angle of the trajectory varied slightly from patient to patient, depending on the frame orientation and location of the frontal burr hole in front of the coronal suture. IC, Internal capsule; SNr, substantia nigra pars reticulata. B, Visualization of all waveforms from 20 sorted units from Patient M (Vim/Vop). Waveforms are 32 samples long and triggered by threshold crossing at sample 8. Individual traces were smoothed and summed into an accumulation buffer. Lighter colors represent a higher density of voltage traces. Traces are labeled by channel followed by unit number. Voltage amplitudes are normalized for display and are therefore represented in arbitrary units.
For both targets, once single-unit recordings had been performed for localization, a 32-channel Pt/Ir microwire (35 μm diameter) array (Ad-Tech Medical Instrument) was passed to the appropriate depth via an outer cannula (Patil et al., 2004), where significant activity was noted with the single-unit electrode. After allowing a few minutes for initial recordings to stabilize, the microwire array was slowly advanced through the cannula. Once the number of clearly distinguishable single units was maximized, the microwire array was left in place. At each electrode depth, the patient was instructed to proceed with the voluntary motor task.
Following completion of the multichannel recording sessions, the microwire array was removed and the DBS treatment electrodes were implanted. As a clinical routine, a brain computed tomography scan was performed within 12 h of the surgery procedure, and in no instance was a hemorrhage or other complication noted. Hence, the clinical risks of temporary placement of the 32 channel microwire array were demonstrated to be very low, as previously reported (Patil et al., 2004). Furthermore, in our research, we have used this electrode in many (N = 72) patients over several years with no post-op evidence of hemorrhage.
Electrophysiological recording.
The 32-channel microwire recordings were performed with a Plexon MAP system. Since this study was performed intraoperatively during electrophysiological mapping of the implantation sites, the recordings for each patient consisted of one or more epochs (up to 8, mean = 3.8), between which the electrode depths were altered. For each recording epoch, single units were sorted offline using custom software developed in-house. Extracted spikes consisted of 32 samples each, sampled at 40 kHz, aligned on the crossing of a linear voltage threshold. Sorting was done by projecting these waveforms into a two-dimensional principal component space. Clusters in principal component space were determined visually and selected by use of a lasso tool to define spatial boundaries. If the electrodes were moved (clearly visible as the microdrive corrupted the recording traces) or threshold was otherwise changed during a recording epoch, the record was broken into separate epochs. Recorded neuronal discharges on a given channel that could not be sorted and isolated as single units were classified as belonging to a multiunit, indicating a collection of individual discharges whose identity could not be firmly established. Figure 1B shows the sorted unit waveforms from a single recording session. We estimated the signal-to-noise ratio (SNR) by dividing the variance of the extracted spike samples by an estimate of the noise variance (Bankman et al., 1993), yielding a mean SNR of 4.69 for all sorted units.
Voluntary motor task.
Patients were placed in a supine, semisitting position in front of a computer monitor. A 5DT Systems Data Glove 5 Ultra haptic glove was placed over the hand contralateral to the microelectrode array. This glove was used to measure flexion/extension of the fingers, sampled at 1 kHz. The average flexion/extension signal, which effectively measured opening/closing of the hand, was used to control the one-dimensional position of the cursor on a video screen, placed directly in front of the patient for high visibility. Patients were trained to modulate the opening and closing of their hand to acquire targets by moving the cursor into a box placed randomly along a horizontal line (Fig. 2A). The required target hold time was 200 ms. Once the target was acquired, the box disappeared for 300 ms before reappearing in a new random position, chosen from a uniform distribution representing the horizontal extent of the screen. Therefore, movements did not strictly alternate between left and right; successive jumps would frequently occur in the same direction.
A, Diagram of the bidirectional hand task. The patient uses graded opening/closing of the hand to actuate a one-dimensional cursor toward randomly appearing targets. B, Example off-line prediction of hand/cursor position. Prediction was obtained using a linear Kalman filter with a 500 s training period, from Patient M (Vim/Vop). Forty-eight units (16 single units, 32 multiunits) were included in the ensemble. Position is given in normalized units.
During the preliminary training/calibration phase, the cursor gain, offset, and target box size were systematically calibrated by the experimenters to compensate for variations in physical ability. Specifically, in patients exhibiting limited hand mobility, the gain from hand movement to cursor motion was increased. In patients with pronounced hand tremor, the size of the target was increased. Finally, the offset was set so that the resting position of the patient's hand corresponded to the middle of the screen. During the recordings that followed, the length of individual motor task sessions varied depending on electrophysiological recording quality and the level of patient fatigue. Figure 2B shows a representative snapshot of the motor task performed by an ET patient. Note the slight 5 Hz tremor that occurred during target hold periods.
Neuronal tuning to target and movement.
All subsequent data analyses were performed using Matlab (MathWorks). We use the term “tuning” to refer to modulation of neuronal firing rates that is correlated to an external parameter and “tuning strength” to refer to the extent of those modulations.
For all sorted single units and multiunits, perievent time histograms (PETHs) of neuronal activity (Awiszus, 1997) were generated using one of two event triggers: (1) the appearance of a new target or (2) movement time. PETHs triggered on target appearance were constructed using a window beginning 0.5 s before each event trigger and ending 1.5 s after, whereas a symmetric 2 s window was used for PETHs triggered on movement time. Movement time was defined as the moment at which the cursor crossed the midpoint between initial cursor position (at target appearance) and the endpoint target position. We chose this standard as a robust definition of movement execution in light of patient tremor and occasional incorrect movements. Regardless of reaction time, this event trigger was locked to movement, being in close proximity to the point of maximum hand velocity before target acquisition. Trials with anomalous movement times <200 ms (premature movement) or >1000 ms (inattention) were discarded. Only neurons with at least 50 valid target acquisition trials were chosen for further analysis.
Neuronal tuning to either target or movement was determined by quantifying the deviation of each PETH from a bootstrap distribution of PETHs generated by uncorrelated triggers. We calculated significance using the one-sample Kuiper's test (Kuiper, 1962; Batschelet, 1981; Zar, 1999), a nonparametric test related to the Kolmogorov–Smirnov (K-S) test (Zar, 1999) but better suited for nonbiased PETH analysis. Unlike the K-S test, Kuiper's test is equally sensitive throughout the distribution, a useful property in scenarios in which the locations of the peak modulations are not known a priori. Variations of the K-S test have been used previously in the significance evaluation of neuronal PETHs (Ghazanfar et al., 2001; Wiest et al., 2005; Gutierrez et al., 2006). In this study, we used Kuiper's test to distinguish an observed distribution of event-triggered spike times from the null hypothesis (uniform probability distribution). Kuiper's test requires the calculation of the maximum positive and negative deviations of the observed PETH cumulative distribution function (CDF) from a uniform distribution CDF (ramp function); the sum of these two deviations is the statistic V: V = max[CDFsample − CDFuniform] + max[CDFuniform − CDFsample].
The Kuiper statistic, K, is a normalized version of V, taking into account the size of the observed sample size N, in this case, the number of binned spikes: (K = VN1/2 + 0.155 + 0.24N−1/2). To distinguish the test statistic Kobs from the null hypothesis, we generated a bootstrapped distribution of 1000 simulated Kuiper statistics (Ksim). Preliminary analysis determined shuffling of spike timestamps to be a suboptimal control; the process eliminates spike autocorrelations from the bootstrap distribution, thereby potentially biasing the evaluation of the observed distribution in favor of significance. Instead, each value of Ksim was calculated using a PETH constructed from the original spike timestamps but processed using a distribution of randomized event triggers; the triggers were drawn uniformly from the time span of the recording session. In other words, in each bootstrapped trial, the timestamps of the events (target appearance, movement) are randomly and independently assigned to decorrelate them from the spiking data. For each sorted unit, the resulting bootstrapped distribution of Ksim was used to produce a p value: p = (number of trials for which Ksimn > Kobs)/(N + 1).
Units were deemed to be tuned to task events (target or movement onset) using the threshold p < 0.05. These units exhibited temporal modulations in firing rate relative to newly appearing targets and/or target-directed movements. Tuning strength was defined as the z-score of the observed PETH relative to the bootstrap distribution.
Directional tuning.
The directional tuning of each sorted unit or multiunit was defined as the difference in neuronal response for leftward versus rightward movements. Significance of directional tuning was determined using the two-sample Kuiper's test. This test applied the calculation of the maximum positive and negative deviations between the spike CDFs for leftward and rightward movements: V = max[CDFsample 1 − CDFsample 2] + max[CDFsample 2 − CDFsample 1].
The Kuiper statistic K was calculated from V using the same equation as the one-sample Kuiper's test, but in the case of the two-sample Kuiper's test, the effective sample size (Neff) replaced N to account for the combined contributions of the individual sample sizes N1 and N2: Neff = (N1N2)/(N1 + N2).
Two PETHs triggered on movement time were generated, one for leftward movements and one for rightward movements. As with the one-sample Kuiper's test, each unit's Kobs statistic was compared with a bootstrapped distribution of Ksim generated from randomized trigger times. Units were deemed to be directionally tuned using the threshold p < 0.05. Tuning strength was defined as the z-score of the observed PETHs relative to the bootstrap distribution.
Tremor sensitivity.
Although concurrent surface EMG recordings were not permitted under our approved experimental protocol, haptic/position tracking have been used repeatedly to analyze tremor (Bardorfer et al., 2001; Su et al., 2003; Vinjamuri et al., 2009), as has accelerometry (Ghika et al., 1993; Grimaldi et al., 2007; Birdno et al., 2008).
For all sorted single units, perievent phase histograms (PEPHs) of neuronal activity were generated using phase of the patient's tremor as an event trigger, similar to the approach of Lebedev et al. (1994). Tremor was determined from hand velocity, and tremor periods within the range of 100–2000 ms (0.5–10 Hz) were analyzed. To exclude the impact of voluntary movements, hand velocity peaks occurring within 250 ms of a movement trigger were excluded. Each tremor period was defined in units of phase, with neuronal spike activity captured into 100 bins of equivalent phase aperture (3.6° each). The zero phase for each cycle was defined by a local maximum in hand velocity. Only sorted units with at least 500 valid tremor periods were chosen for further analysis.
Each resulting PEPH was a measurement of neuronal firing rate with respect to tremor phase. The one-sample Kuiper's test, in addition to possessing uniform sensitivity, is also rotationally invariant, meaning that the arbitrary choice of zero phase has no effect on the assessment of statistical significance. For analysis of tremor tuning, we generated a bootstrapped distribution of 1000 simulated Kuiper statistics (Ksim); each was calculated using a PEPH constructed from trials whose binned spike counts were circularly rotated by uniformly random phase offsets. For each analyzed unit, the bootstrapped distribution was used to produce a p value. Units were deemed to be tremor associated (tuned) using the threshold p < 0.05. Tuning strength was defined as the z-score of the observed PEPH relative to the bootstrap distribution.
Oscillatory neurons.
For all sorted units with at least 1000 extracted spikes, we used Welch's method (Oppenheim and Schafer, 1975) with eight nonoverlapping segments to determine the spike train autopower spectral density. The power spectra were smoothed using a 0.5 Hz rectangular sliding window. For each unit, the peak autopower frequency was determined in the 1–25 Hz range, with frequency content <1 Hz discarded for the remainder of the analysis. For the peak frequency, we determined the SNR by dividing peak power by the mean power (assessed from 1 Hz up to the Nyquist frequency of 500 Hz). For the purpose of comparison, the same spectral analysis was performed on hand acceleration traces for all recorded sessions.
Preliminary analysis indicated a functional separation of peak frequencies at ∼2.5 Hz. Units with a peak power frequency <2.5 Hz tended to be dominated by low-frequency power and were therefore judged not to be sufficiently oscillatory in a physiologically relevant frequency range. As in previous studies (Lenz et al., 1988; Amtage et al., 2008), only units with a peak SNR >2 were classified as oscillatory. Sharpness of the peaks in either the spike train or hand acceleration autopower spectra was determined by calculating the maximum power concentrated in a 1 Hz band within the physiologically relevant 2.5–7.5 Hz window. “Peakedness” was defined as the ratio of the power in this band relative to the total power in the 1–25 Hz band.
In addition to established linear methods, we developed a heterodyne method for detecting nonstationary or wide-bandwidth synchronized activity between each neuron's firing rate and associated hand velocity. This method, inspired by the frequency shifting scheme used in radio frequency transceivers, was applied to all sorted units that fulfilled the selection criteria for both target tuning analysis and oscillatory analysis. The spike train and hand velocity recordings were first bandpass filtered (fourth-order Butterworth, zero phase method) between 2 and 12 Hz, leaving both signals with negligible 0 Hz (DC) energy. The two signals were then multiplied, yielding a third time series from which to extract spectral energy. Because multiplication in the time domain is equivalent to convolution in the frequency domain (and vice versa), any synchronous frequency-modulated components in both neuronal firing rate and hand velocity are transferred to DC. The ratio of spectral energy (Eobs) from 0 to 0.125 Hz (signal) over that from 0.25–2 Hz (baseline) was considered as a metric of heterodyne correlation between neuronal activity and hand movement. To provide a control for this estimate, a bootstrap distribution (Esim) was generated by shuffling spike timestamps 1000 times and repeating the above analysis. Since all low-frequency information is filtered out before multiplication, shuffling timestamps was determined to be a suitable control. For all single units, the resulting bootstrapped distribution was used to produce a p value. Units were deemed to be heterodyne-tuned to tremor if the heterodyne correlation between neuronal activity and movement was significant using the threshold p < 0.05. Heterodyne tremor tuning strength was defined as the z-score of the observed spectra relative to the bootstrap distribution.
Efficacy of neuronal recordings for kinematic predictions.
Prediction algorithms from the BMI literature were applied to subcortical neuronal populations to extract behavioral parameters. Several algorithms were tested, including the linear Kalman filter, unscented Kalman filter, and the Wiener filter. Figure 2B shows an example off-line prediction for a 30 s window of task performance.
Since all three algorithms achieved the same approximate fidelity in preliminary testing, we chose the Wiener filter for further analysis due to its computational simplicity and extensive presence in the BMI literature (Wessberg et al., 2000; Carmena et al., 2003; Patil et al., 2004). Individual Wiener filters were fit by binning neuronal data into 100 ms time slices with 10 causal lags and regressing against recorded hand position. Model training was performed by the random selection of 50% of these time slices; predictions were then made on a distinct random 25%. This process was repeated with 100 draws of fit and predict time slices. Correlation coefficient (R) between predicted hand position and actual hand position was measured for each of the draws; the mean correlation coefficient was reported for each recording session. Offline prediction results are reported for all sessions with at least 50 presented targets.
We generated neuron dropping curves for selected sessions (Wessberg et al., 2000) by drawing random subsets from the neuronal ensemble. For each subset ensemble size N, we performed 1000 draws of random ensemble subset and Wiener filter fit and prediction; the R values for these draws were averaged to form a smooth neuron dropping curve. Following Wessberg et al. (2000), the resulting curve was then fit to the following hyperbolic function to extrapolate the performance results to larger ensemble sizes: R2 = cN/(1 + cN).
Neuronal synchrony.
The use of simultaneous ensemble recordings allows for the analysis of pairwise synchrony between neurons. To determine the statistical significance of the synchrony between two neurons, we analyzed the cross-correlation peak between pairs of spike trains. Pairs were analyzed if they each contained at least 100 spikes and corresponded to a session with at least 50 targets. The cross-correlation coefficient was first calculated for the observed spike trains of the two neurons, then smoothed using a 5 ms rectangular sliding window. The observed test statistic Cobs was defined as the peak coefficient in the ±10 ms time lag range. Bootstrap simulations (n = 1000) of the two spike trains were generated by convolving the spike trains with a Gaussian kernel (σ = 250 ms) and then generating new spike trains via an inhomogeneous Poisson process. The smoothing filter was used to extinguish correlated high-frequency content in the bootstrap distribution while maintaining low-frequency correlation in mean firing rate. Each of these bootstrap simulations was used to produce a cross-correlation coefficient Csim. The bootstrapped distribution was used to produce a p value, and neuron pairs were deemed to be significantly synchronous using the threshold p < 0.05.
To visualize the time dependence of pairwise synchrony, we generated joint peristimulus time histograms (JPSTHs) for neuron pairs, as originally proposed by Aertsen et al. (1989). Our JPSTHs were adjusted by subtracting the shift predictor histogram and normalizing (bin-by-bin) by the standard deviation, a procedure referred to by Aertsen et al. (1989) as the “true normalization” of the JPSTH.
Results
A total of 25 DBS implantation patients were examined. In these patients, we simultaneously recorded from ensembles of up to 23 well isolated neurons from either Vim/Vop or STN, depending on the site of electrode location. Recording sessions varied substantially in terms of duration and target acquisition rate, as limited by individual patient pathology and motivation. Neurons from these subcortical areas were classified by oscillatory firing patterns and tuning to target, movement, direction, and tremor. Moreover, neuronal ensemble data served as input for an offline linear prediction model to reconstruct cursor position. Finally, neuronal pairs were analyzed for evidence of functional synchrony.
STN cells (N = 168) exhibited a higher (p < 0.01, Mann–Whitney U test) mean firing rate than Vim/Vop cells (N = 83): 15.8 ± 1.95 Hz and 11.7 ± 1.02 Hz, respectively (mean ± 1 SE in both cases). In both subcortical areas, we found substantial populations of oscillatory neurons, as well as neurons strongly tuned to target, movement, direction, and tremor. Furthermore, neurons in both subcortical areas tended to show tuning to multiple parameters (Table 1) rather than belonging to disjoint sets. For example, the number of Vim/Vop cells tuned to both target and tremor was higher than would be expected under statistical independence (two-tailed Fisher's exact test, p < 0.05). At the ensemble level, a substantial number of analyzed cell pairs were found to exhibit synchrony. Curiously, all tremor-associated neurons exhibited synchrony within the recorded neuronal ensemble.
Pairwise classifications for single units
Neuronal tuning to target and movement
Both Vim/Vop and STN neurons represented target appearance and movement onset (Table 2). Of all single units tested, 29.2% of 168 Vim/Vop cells and 22.9% of 83 STN cells were found to be tuned to target appearance. Both of these percentages represent significant populations (binomial test, p ≪ 0.001 in both cases). Figure 3A shows example PETHs for three highly responsive neurons.
Behavioral tuning of subcortical neurons
Example PETHs. A, Strongly tuned units to target appearance. i, Vim/Vop cell, Patient M, 465 trials. ii, Vim/Vop cell, Patient M, 375 trials. iii, STN cell, Patient H, 310 trials. B, Strongly tuned units to movement time. i, Vim/Vop cell, Patient M, 390 trials. ii, Vim/Vop cell, Patient M, 310 trials. iii, STN cell, Patient H, 201 trials. C, Strongly direction tuned units. i, STN cell, Patient H, 201 trials. ii, Vim/Vop cell, Patient M, 416 trials. iii, Vim/Vop cell, Patient M, 394 trials. For A, reported time is relative to target appearance. For B and C, reported time is relative to movement time. All plots were smoothed using a Gaussian kernel, σ = 40 ms.
As explained above, trials with reaction times outside the 200–1000 ms range were discarded for movement tuning, and only sessions with at least 50 valid trials were subjected to further statistical analysis. Because of the additional reaction-time criterion, fewer single units were analyzed for movement tuning than for target tuning. Of these, 34.7% of 75 Vim/Vop cells and 42.3% of 26 STN cells were found to be tuned to movement. Both of these percentages represent statistically significant populations (binomial test, p ≪ 0.001 in both cases). Figure 3B shows example PETHs for three highly responsive neurons.
We also found a strong positive correlation between the strength of target appearance tuning and that of movement tuning, for both Vim/Vop and STN cells. When controlling for the number of session trials, target tuning strength significantly predicted movement tuning strength (β = 0.70, p ≪ 0.001 for Vim/Vop; β = 0.71, p ≪ 0.001 for STN). This result is consistent with Table 1, which indicates that a larger-than-expected number of neurons in both subcortical areas were tuned to both target and movement.
A portion of the correlation between target tuning and movement tuning may be explained by a tight temporal offset between target appearance and movement time. However, visual inspection of some tuned units indicated a clear decoupling of the neural encoding of target appearance and movement. Sorted raster plots from example neurons are shown in Figure 4, A and C. From these, we derived the color maps in Figure 4, B and D, each showing two clear bands of increased spike density. In both panels, the vertical bands are independent of movement time and are clearly related to target appearance (∼450 ms postappearance). The diagonal bands have a near-unity slope, indicating a clear time-locked relationship between neuronal activity and movement time. For both units, the second peak in firing rate occurred ∼300 ms after the defined movement time. From these data, it can be concluded that these neurons were tuned to both target appearance and movement; they modulated their firing rates in relation to both events.
Separation between single-unit response to target appearance and movement time. Time along the x-axis is relative to target appearance. A and C show spike raster plots relative to target appearance, with individual trials sorted by movement time (red circles). B and D show smoothed color plots for the same two units as A and C, using the relative timing of all spikes and the movement time of their corresponding trials. The dashed line (unity slope) depicts movement time. For generation of the color plots, the data were smoothed using a two-dimensional Gaussian kernel, with σ = 20 ms along the x-axis and σ = 40 ms along the y-axis. Both units are Vim/Vop cells from Patient M.
Modest differences were seen in the aggregate response patterns of Vim/Vop and STN cells classified as responsive to either target or movement (Fig. 5). Both cell types exhibited a mean response that peaked following target appearance (Fig. 5A); STN cells peaked later on average. The mean response of the Vim/Vop cells peaked immediately before movement while that of the STN cells peaked concurrently with movement (Fig. 5B). Differences in the relative lags for Vim/Vop and STN neuronal activation likely reflect the position of thalamic and STN neurons in the network hierarchy of motor control (Marsden et al., 2001; Guillery and Sherman, 2002; Gradinaru et al., 2009). The motor regions of the thalamus are more involved with intention, with signals arriving before motor cortex activation, whereas the collaterals from motor cortex to STN deliver signals at the time of motor activation.
A, B, Mean normalized PETHs for all responsive units (p < 0.05) from both subcortical areas using two event triggers: target appearance (A), and movement time (B). Reported time is relative to the event trigger. Before aggregation, individual PETHs were convolved with a Gaussian kernel (σ = 40 ms) and normalized relative to mean firing rate. For both panels, the differences between the mean Vim/Vop and STN responses were statistically significant (χ2 test, p ≪ 0.001).
For both aggregates, the differences between the mean Vim/Vop and STN responses were statistically significant (χ2 test, p ≪ 0.001). However, similar proportions of Vim/Vop and STN cells were tuned to target appearance; the same was also true for movement tuning (two-tailed Fisher's exact test, p > 0.05 in both cases).
Directional tuning
Another metric of interest for the behavioral responsiveness of subcortical neurons was directional tuning; Table 2 gives the directional tuning results for both single units and multiunits. Of all tested single units, 25.3% of 75 Vim/Vop cells and 19.2% of 26 STN cells were found to exhibit directional tuning. Both of these percentages represent statistically significant populations (binomial test, p < 0.001 for Vim/Vop, p < 0.01 for STN). Figure 3C shows example PETHs for three strongly tuned neurons. Similar proportions of Vim/Vop and STN cells were tuned to direction (two-tailed Fisher's exact test, p > 0.05).
Despite the clear separation in the neuronal response to leftward and rightward movements, note the transient regions of convergence that occurred in Figure 3C. In Figure 3Ciii, for example, the neuronal responses to each direction converged just before movement. For many tuned neurons in both Vim/Vop and STN, the degree of directional modulation varied throughout the temporal window.
For both Vim/Vop and STN cells, we found strong positive correlations between directional tuning strength and the strength of both target tuning and movement tuning. When controlling for the number of session trials, target tuning strength significantly predicted directional tuning strength (β = 0.16, p < 0.05 for Vim/Vop; β = 0.38, p < 0.01 for STN). Similarly, movement tuning strength significantly predicted directional tuning strength (β = 0.24, p < 0.01 for Vim/Vop; β = 0.43, p < 0.05 for STN). The latter finding is consistent with the Vim/Vop pairwise classification result in Table 1.
Properties of multiunits
Whereas single units are identifiable as distinct neurons, a multiunit is likely comprised of distant neurons with lower SNR spike profiles. From Table 2, it can be seen that a significant population of analyzed multiunits were tuned to target, movement, and direction (binomial test, p < 0.01 in all cases). A substantial number of these tuned multiunits were found on the same recorded channel as tuned sorted units. Furthermore, when controlling for the number of session trials, the target tuning strength of single units significantly predicted the target tuning strength of same-channel multiunits (β = 0.18, p < 0.01). This confirms the presence of correlated tuning in nearby neurons. Thus, a substantial amount of encoded information was present in subcortical multiunits, arguing for the potential inclusion of these signals in future analyses of ensemble activity.
Tremor sensitivity
To identify potentially pathological neurons within the recorded subcortical populations, we analyzed the tremor sensitivity of single units using the discussed PEPH approach; the results are given in Table 3. Of all single units tested, 12.4% of 169 Vim/Vop cells and 15.9% of 82 STN cells were found to be correlated to observable hand tremor. Both of these percentages represent statistically significant populations (binomial test, p < 0.001 in both cases). Figure 6 shows example PEPHs for three strongly tremor-sensitive neurons. These results demonstrate that for highly tuned units, the dependence of spike rate on tremor phase remained stable throughout the recording session (Fig. 6), even if the mean firing rate varied substantially Figure 6, A and C. Similar proportions of Vim/Vop and STN cells were tuned to tremor (two-tailed Fisher's exact test, p > 0.05).
Tremor tuning of subcortical neurons
Example PEPHs triggered on hand tremor phase for strongly tremor tuned units. A, Vim/Vop cell, Patient M, 4478 tremor periods. B, Vim/Vop cell, Patient M, 1259 tremor periods. C, Vim/Vop cell, Patient M, 4642 tremor periods. Zero phase is aligned to local peaks in hand velocity. All panels show individual PEPHs for each temporal quarter of their respective sessions, overlaid with the aggregate PEPH for the session as a whole (ALL). All plots in this figure are smoothed using a Gaussian kernel, σ = 15°.
For Vim/Vop cells (but not STN cells), we found a positive correlation between the strength of tremor tuning and that of directional tuning. When controlling for the number of session trials, directional tuning strength significantly predicted tremor tuning (β = 0.34, p < 0.05). However, we found no relationship between tremor tuning and either undirected target or movement tuning (p > 0.05 for all cases).
However, these results do not distinguish whether these tremor-tuned neurons are involved in a pathological mechanism that causes tremor or merely reflect somatosensory signals indicative of tremor.
Oscillatory behavior
To explore neuronal oscillations in Vim/Vop and STN and their relationship to patient pathology, we inspected the autopower spectra of single-unit spike trains for strong frequency peaks, yielding peak frequency and SNR (Fig. 7A). Of all tested single units with SNR > 2, the distribution of peak frequencies showed a clear bimodal distribution with a border between low- and high-frequency oscillations at ∼2.5 Hz (Fig. 7B).
A, Distribution of peak frequencies for the spike train autopower spectra of 397 analyzed single units. Units with peak frequency >10 Hz (17.4% of all units) are not shown in this figure; they are distributed with near uniformity in the 10–25 Hz range with SNR mostly below the classification threshold of 2. B, Histogram showing the number of oscillatory (SNR > 2) units at each peak frequency. A clear separation between two subpopulations is clear at ∼2.5 Hz.
Spike train spectra from Vim/Vop and STN neurons also tended to possess large amounts of energy at low frequencies, suggestive of 1/f (pink) noise. This power-law distribution has been described for cortical neurons as a stochastic process (Davidsen and Schuster, 2002), but may also be related to slow modulations of patient attention and arousal. The observed distribution may explain the 2.5 Hz trough (Fig. 7) that separates neurons dominated by 1/f noise from those exhibiting strong oscillations in the tremor-relevant frequency range (2.5–7.5 Hz). Only neurons with sufficient power within this frequency range were eligible to be classified as oscillatory.
The oscillatory classification results are shown in Table 4; 21.5% of 274 Vim/Vop cells and 17.9% of 123 STN cells were classified as oscillatory. No difference was seen in the proportions of oscillatory Vim/Vop and STN cells (two-tailed Fisher's exact test, p > 0.05). Figure 8 shows example interspike interval (ISI) plots for three highly oscillatory cells. Note that all three ISI histograms exhibit some degree of bimodality, indicative of periodic bursting behavior.
Oscillatory Units
Example ISI histograms for highly oscillatory neurons. A, Vim/Vop cell, Patient W. B, STN cell, Patient V. C, Vim/Vop cell, Patient J.
Figure 9A shows the smoothed autopower spectra of spike trains for all analyzed single units, with each individually normalized horizontal trace corresponding to a distinct unit. From this figure, one can visually identify some of the highly oscillatory units as well as observe the congruity between multiple units from the same patient. The difference between the mean normalized spectra for Vim/Vop and STN cells (Fig. 9B) is statistically significant (χ2 test, p ≪ 0.001). From Figure 9B, it is clear that STN cells tended to concentrate power at a lower frequency (3 rather than 4 Hz).
A, Smoothed autopower of spike trains for all sorted units with sufficient spike count. The autopower spectra (determined by Welch's method) of each horizontal trace, corresponding to a distinct unit, has been individually smoothed (using a 0.5 Hz sliding window) and normalized. The horizontal traces are grouped by patient and recording area. The color bar represents arbitrary units of normalized energy density. B, Mean autopower spectra from all analyzed Vim/Vop and STN cells. Spectra are individually smoothed using a 0.5 Hz sliding window and normalized before aggregation. Each of the dashed traces represents a bootstrap simulation in which only half of all analyzed units are aggregated. The difference between the mean Vim/Vop and STN spectra is statistically significant (χ2 test, p ≪ 0.001).
The pairwise classification results in Table 1 reject the notion that oscillatory neurons and behaviorally tuned neurons form disjoint sets. Furthermore, we found no relationship between spike autopower peakedness and the strength of any of the three (target, movement, direction) behavioral tuning metrics (p > 0.05 for all cases, for both Vim/Vop and STN). The lack of a clear anticorrelation suggests that the sets of behavioral neurons and oscillatory neurons are far from disjointed. Instead, they appear to exist as overlapping populations.
Our next analysis intended to uncover a relationship between strong oscillatory neuronal patterns and observable hand tremor. However, we did not find any clear relationship. Linear regression analysis revealed no relationship between peak frequency (2.5–7.5 Hz range) of spike train autopower spectra and corresponding hand acceleration autopower spectra (p > 0.05 for both Vim/Vop and STN). No relationship was found between the sharpness of the two spectra for Vim/Vop neurons (p > 0.05), but we did observe a marginally significant positive relationship for STN neurons (β = 0.72, p = 0.038). Figure 10 shows overlaid spectra for the spike train autopower and hand acceleration autopower of three highly oscillatory units. For all three cells (representative of the population as a whole), the peak frequencies do not coincide. However, we did find a marginally significant (β = 0.24, p = 0.055) correlation between spike autopower peakedness and tremor tuning strength for Vim/Vop cells (p > 0.1 for STN cells). These findings call into question the presumed causal linear relationship between the two, suggesting the possibility of an elusive nonlinear relationship.
Comparison of spike train autopower and hand acceleration autopower for three highly oscillatory units. A, Vim/Vop cell, Patient W. B, STN cell, Patient V. C, Vim/Vop cell, Patient J. Note that for all three cells, peak frequency does not coincide.
Our heterodyne decoding analysis further explored this relationship by applying a nonlinear frequency shifting approach to the autopower spectra. Of 239 analyzed single units, 13.3% of Vim/Vop cells and 17.3% of STN cells were found to be tremor associated via heterodyne decoding (Table 5). Both of these percentages represent statistically significant populations (binomial test, p ≪ 0.001), but the difference between them is not significant (two-tailed Fisher's exact test, p > 0.05).
Heterodyne tremor tuning
Heterodyne decoding is a more sensitive method for detecting tremor correlations than spectral peak analysis if the tremor signal is wide-bandwidth or prone to phase changes. Indeed, this schema may better serve to explain the relationship between the oscillatory activity of neurons and observed tremor. In fact, this can explain the similarity in the proportions of tuned neurons in Tables 3 and 5. Furthermore, Vim/Vop firing indicated a strong positive correlation between tremor tuning strength, identified using PEPHs, and heterodyne tremor tuning strength (β = 0.31, p < 0.001). This relationship was marginally significant in STN cells (β = 0.28, p = 0.086). These findings are consistent with the pairwise classification results in Table 1, which indicated a higher than expected joint classification for the two tremor tuning analyses for Vim/Vop cells.
Efficacy of neuronal recordings for kinematic predictions
We also performed off-line predictions of cursor motion using the recorded ensembles. The correlation coefficient (mean ± 1.98 SE) for each of the sessions is shown in Figure 11. Although the predictions varied greatly across sessions and patients, the results compared favorably with our previous study (Patil et al., 2004). The best session for each subcortical area (Vim/Vop, STN) was chosen for further analysis, and neuron dropping curves were generated for these two sessions and fitted to a hyperbolic function (Wessberg et al., 2000). Extrapolation of the hyperbolic fit produced estimates of the approximate ensemble sizes required to achieve R2 = 0.9: 106 Vim/Vop neurons or 397 STN neurons.
Dependence of off-line BMI predictions on neuron ensemble size. Each data point corresponds to a recording session. Correlation coefficient (R) is indicated as mean ± 1.98 SE. The best session for each subcortical area (Vim/Vop, STN) was chosen for further analysis. Neuron dropping curves are shown for each of these sessions and fitted to a hyperbolic function (Wessberg et al., 2000). Extrapolation calculations are presented in Results.
Neuronal synchrony
We analyzed neuronal synchrony in pairs of sorted units and investigated how its prevalence varied across subcortical areas; the results are given in Table 6. Using the cross-correlation approach, 43.0% of Vim/Vop pairs and 25.8% of STN pairs were found to be significantly synchronous. Both of these percentages represent significant populations (binomial test, p ≪ 0.001 for both cases). Figure 12A shows example cross-correlation plots for three highly synchronous pairs, while Figure 12B shows the normalized JPSTH for the same three pairs. Whereas Figure 12B, i and ii, clearly shows temporal synchronization along the diagonal (and off-diagonals), the same result is not visually discernible in Figure 12Biii.
Pairwise neuronal synchrony
A, Example plots of cross-correlation coefficient for three highly synchronous neuron pairs. Bootstrapped simulations of coefficient shown in red. B, Example normalized JPSTHs for the same three neuron pairs. i, Pair of Vim/Vop cells, Patient I. ii, Pair of STN cells, Patient H. iii, Pair of Vim/Vop cells, Patient Q. For B, reported time is relative to target appearance.
We found highly significant differences between the Vim/Vop and STN in terms of the proportions of synchronous pairs. A significantly higher proportion of Vim/Vop pairs were synchronous than STN pairs (two-tailed Fisher's exact test, p ≪ 0.001).
It has been reported that the level of tremor in parkinsonian patients is positively correlated to the degree of pairwise synchrony among STN cells (Levy et al., 2000). To test the relationship between tremor tuning and local synchrony, we compared the subpopulation of both Vim/Vop and STN neurons that were synchronous with at least one other neuron in their respective ensembles to the subpopulation of neurons tuned to hand tremor (PEPH method). Only neurons fulfilling the criteria of both individual analyses were considered. The results are shown in Table 7. The observed proportions are significantly different (two-tailed Fisher's exact test, p < 0.01), indicating a clear interaction between tremor tuning and local synchrony. Only units synchronized to at least one other unit were tuned to tremor, whereas no unsynchronized units were tuned to tremor.
Comparison of synchrony (cross-correlation method) and tremor tuning (PEPH method)
Discussion
In this study, we analyzed ensemble activity of human subcortical neurons (either Vim/Vop or STN) in 25 patients, recorded in patients who performed visually guided hand movements. To our knowledge, this constitutes the largest sample of human subcortical ensemble recordings to date. We quantitatively evaluated the representation of motor parameters in these neurons, as well as activity thought to be related to pathological states—tremor sensitivity, oscillations, and pairwise synchrony. The present work also supports our previous proposition that a sufficiently large ensemble of subcortical neurons could enable a motor BMI (Patil et al., 2004). We further propose that chronic subcortical microelectrode technology could serve as the basis of a new generation of neuroprosthetic devices aimed at both monitoring and actively reducing oscillatory firing and elevated levels of network synchrony.
Subcortical encoding of behavior
Neurons in both subcortical areas (Vim/Vop and STN) were found to encode features of patient motor behavior, both voluntary (target tracking) and involuntary (tremor). During voluntary behavior, a substantial number of neurons were found to be tuned to target appearance, movement onset, and movement direction (Table 2). This finding likely indicates that neurons in both structures are broadly tuned across multiple modalities and muscle groups. This observation is consistent with previous studies in which STN neurons have been reported to have large receptive fields that respond to multiple joints (Abosch et al., 2002), with 40% of STN cells tuned to movement in a simple two-dimensional joystick task (Williams et al., 2005). It has been reported that 42% of STN neurons respond to passive movement of either arm or leg, and of these, 25% responded to multiple joint movements (Theodosopoulos et al., 2003), whereas 51% of Vim/Vop neurons are tuned to sensory stimuli, with an overlapping 10% of these cells tuned to volitional movements (Lenz et al., 1990).
Some thalamic and STN neurons exhibited bimodal encoding of target and movement (Fig. 4), demonstrating what appears to be superposition of two independent neural representations. Furthermore, we observed a significant overlap between tuning to target appearance and to movement onset in both structures. Additionally, directional tuning strength was positively correlated to movement tuning strength, indicating that individual neurons encoded both parameters by exhibiting both common-mode and differential activity relative to movement onset (Fig. 3). Thus, neurons from both subcortical areas exhibited rate modulations based on a broad superposition of task parameters (Table 1).
Our analysis of involuntary motor activity (tremor) yielded phase histograms demonstrating clear phase-locking between neuronal activity and tremor periods (Fig. 6). Overall, we report significant populations of tremor tuned cells (12.4% for Vim/Vop, 15.9% for STN). The literature reports a large range in the prevalence of tremor-related cells. For STN, researchers have reported 11% (Magariños-Ascone et al., 2000), 19% (Rodriguez-Oroz et al., 2001), and 52% (Amtage et al., 2008). For Vim/Vop, researchers have reported 34% (Lenz et al., 1988), 35.6% (Zirh et al., 1998), and 51% (Hua and Lenz, 2005). These disparities are probably due to differences in recording parameters and classification methodology. For example, whereas many researchers have classified tremor tuning using linear coherence between spike train and recorded EMGs, we used kinematic hand position recorded by a haptic glove.
We found substantial overlap between tremor-tuned neurons and those tuned to parameters of voluntary motor behavior (Table 1). Moreover, we also observed a positive correlation between tremor tuning strength and directional tuning strength (Vim/Vop cells only). These results are consistent with earlier studies (Magariños-Ascone et al., 2000; Rodriguez-Oroz et al., 2001) that reported a large proportion of tremor-related STN cells to be simultaneously related to voluntary movements either through motor or sensory loops.
Elusive relationship between oscillatory activity and tremor
Within a tremor-relevant frequency range (2.5–7.5 Hz), we observed a substantial population of Vim/Vop and STN cells exhibiting strong oscillations (Table 4). Furthermore, we found that the mean autopower spectra for recorded Vim/Vop cells had a higher peak frequency than that of recorded STN cells (Fig. 9B). This finding is consistent with reports of higher mean tremor frequency in ET patients (4–12 Hz) than in PD patients (3–6 Hz) (Deuschl et al., 1998).
Whereas Rodriguez-Oroz et al. (2001) reported that oscillatory cells in STN did not represent movements, we found moderate overlap (Table 1) between Vim/Vop and STN neurons exhibiting oscillations and those tuned to voluntary behavioral parameters (target, movement, direction). This important finding indicates that there is not a strict dichotomy between pathological neurons and those encoding motor signals.
Similarly, we found overlap but no significant correlation between the presence of oscillatory patterns and tremor tuning in both subcortical areas. This finding corroborates the claims of earlier studies of ventral thalamus and STN (Magnin et al., 2000; Rodriguez-Oroz et al., 2001), in which the sets of tremor-related neurons and oscillatory neurons show modest intersection. In other words, not all tremor-related neurons exhibited oscillations, and some oscillatory neurons exhibited no clear association with tremor.
Furthermore, spectral peak detection methods revealed no relationship between spike train and hand acceleration (Fig. 10), though there may be an elusive nonlinear or nonstationary relationship between neuronal oscillations and hand tremor. We hypothesized that for some cells, spike train and tremor comodulated with relatively large bandwidth. Our heterodyne decoding results suggest that a substantial number of cells in both subcortical areas may be tremor-tuned in this manner (Table 5). Overall, no consensus in the field has been reached regarding the definitive relationship between oscillatory behavior and tremor tuning, and our results support the idea that pathological oscillations are idiopathic across both neurons and patients.
Network synchrony
Dopamine depletion in PD has been reported to promote neuronal synchrony within the basal ganglia (Heimer et al., 2006). In this study, a high percentage of neuronal pairs from both subcortical areas exhibited synchronous behavior (Table 6). This is generally consistent with reports of pairwise synchrony in the majority of analyzed STN pairs (Levy et al., 2000, 2002). However, we observed high levels of functional synchrony in neuronal pairs much further apart than the submillimeter separation in Levy et al. (2000, 2002). Furthermore, Vim/Vop neuron pairs (from ET patients) exhibited significantly higher synchrony than STN pairs.
Remarkably, only cells exhibiting synchrony with another cell in the ensemble were found to be tremor-tuned; no unsynchronized cells were tremor-tuned (Table 7). This finding corroborates reports that the prevalence of STN pairs synchronized at high frequencies is correlated to the degree of parkinsonian tremor (Levy et al., 2000) and that patients without observable tremor do not exhibit high-frequency STN synchrony (Levy et al., 2002). Others have reported correlations between STN synchrony and bradykinesia/rigidity (Weinberger et al., 2009). Our findings have extended Levy et al.'s (2000, 2002) results to combinatorial pairs in neuronal ensembles from both STN and Vim/Vop. It remains unclear to what degree network synchrony is indicative of tremor pathology, although the effective disruption of unstable network activity by DBS stimulation suggests a connection.
Chronic subcortical ensemble recordings
Acute intraoperative recordings exhibit instability and other limitations, whereas long-term recordings from subcortical structures will offer greater signal quality following the recovery of the electrode-tissue interface. Thousands of subcortical DBS implantations are performed every year with minimal risk (Bronstein et al., 2011). Furthermore, long electrode tracts have been shown to yield more stable extracellular recordings (Porada et al., 2000; Krüger et al., 2010) by mitigating the problem of electrode micromotion seen in cortical implants.
We suggest that chronic subcortical ensemble recordings may bring about the viability of subcortical BMI systems, first discussed in our previous study (Patil et al., 2004). Our reported offline prediction results from our best Vim/Vop and STN sessions are comparable to those reported from selected rhesus macaque cortical regions with similar neuron count (Wessberg et al., 2000; Carmena et al., 2003). However, hyperbolic extrapolation (Wessberg et al., 2000) suggests an ensemble size greater than 100 to achieve prediction fidelity above R2 = 0.9. This would of course necessitate the design of microelectrode arrays for cannula-based implantation with more recording sensors and demonstrable long-term safety and recording efficacy.
We propose that clinical studies using chronic ensemble recordings in humans will permit both the continued study of the neurophysiological mechanisms involved in motor control as well as long-term monitoring of pathological activity. Chronic recordings will facilitate continuous examination of changes in tuning, oscillations, and synchrony as a function of the patient's symptomatic state both on and off treatment. Specifically, this wealth of electrophysiological data may well be used to instruct the improvement of closed-loop DBS systems that are currently in initial stages of development (Rosin et al., 2011; Rouse et al., 2011).
We have shown that within the STN and Vim/Vop thalamus, there is an idiopathic mixture of pathology and behavior tuning; these sites are worthy of more targeted treatment than the current dominant approach of delivering high-frequency DBS through macroelectrodes. Subcortical ensembles remain an untapped resource, with the potential to advance both neuroscience and neurorehabilitation alike.
Footnotes
This work was supported by funding from the National Science Foundation Graduate Research Fellowship award DGE-1106401-004 to A.M.F.; the VA Merit Review Award, and NIH Grants R21NS066115 and RO1AG037599 to D.A.T.; and by Defense Advanced Research Projects Agency Grant N66001-06-C-2019, NIH Grant R01NS073125, and the NIH Director's Pioneer Award DP1OD006798 to M.A.L.N. The content is solely the responsibility of the authors and does not necessarily represent the official views of the Office of the NIH Director or the NIH. We thank Daniel Clayton for surgical assistance and expertise, Susan Halkiotis for manuscript editing and preparation, and Joseph E. O'Doherty for technical, scientific, and programming help.
The authors declare no financial conflicts of interest.
- Correspondence should be addressed to Miguel A. L. Nicolelis, 311 Research Drive, Bryan Research Building, Room 327, Duke University, Durham, NC 27710. nicoleli{at}neuro.duke.edu
