Abstract
Functional neurosurgery has afforded the opportunity to assess interactions between populations of neurons in the human cerebral cortex and basal ganglia in patients with Parkinson's disease (PD). Interactions occur over a wide range of frequencies, and the functional significance of those >30 Hz is particularly unclear. Do they improve movement, and, if so, in what way? We acquired simultaneously magnetoencephalography and direct recordings from the subthalamic nucleus (STN) in 17 PD patients. We examined the effect of synchronous and sequential finger movements and of the dopamine prodrug levodopa on induced power in the contralateral primary motor cortex (M1) and STN and on the coherence between the two structures. We observed discrete peaks in M1 and STN power at 60–90 Hz and at 300–400 Hz. All these power peaks increased with movement and levodopa treatment. Only STN activity at 60–90 Hz was coherent with activity in M1. Directionality analysis showed that STN gamma activity at 60–90 Hz tended to drive gamma activity in M1. The effects of levodopa on both local and distant synchronization at 60–90 Hz correlated with the degree of improvement in bradykinesia-rigidity as did local STN activity at 300–400 Hz. Despite this, there were no effects of movement type, nor interactions between movement type and levodopa in the STN, nor in the coherence between STN and M1. We conclude that synchronization at 60–90 Hz in the basal ganglia cortical network is prokinetic but likely through a modulatory effect rather than any involvement in explicit motor processing.
Introduction
There has been considerable recent interest in the exaggerated beta band (13–35 Hz) oscillations found in the cortico-basal ganglia loop of patients with untreated Parkinson's disease (PD) (Uhlhaas and Singer, 2006; Hammond et al., 2007; Eusebio and Brown, 2009; Weinberger et al., 2009). Mounting evidence suggests that this activity may contribute to slowness and stiffness in PD (Kühn et al., 2006, 2008, 2009; Weinberger et al., 2006; Ray et al., 2008; Chen et al., 2010). In contrast, less is known about oscillatory dynamics at higher frequencies in the cortico-basal ganglia loop. Several rhythmic activities have been reported and are increased after treatment with the dopamine prodrug levodopa and with voluntary movement. Within the gamma frequency band, discrete peaks may occur centered between 60 and 90 Hz in local field potential (LFP) activity and in the coherence between different nodes within the basal ganglia loop (Brown et al., 2001; Cassidy et al., 2002; Williams et al., 2002; Fogelson et al., 2005; Alonso-Frech et al., 2006). An additional discrete peak in the subthalamic nucleus (STN) LFP may be seen in the high gamma (150–350 Hz) band (Foffani et al., 2003; Kane et al., 2009), which shifts in frequency with dopaminergic therapy (Ozkurt et al., 2011). These activities in the gamma and high gamma bands have been considered prokinetic (Brown, 2003; Foffani and Priori, 2006; Ozkurt et al., 2011). However, only some features of (high) gamma band activity represent a prokinetic neurophysiological marker of levodopa-induced motor improvement in PD (López-Azcárate et al., 2010; Ozkurt et al., 2011). Our first aim was therefore to test the hypothesis that gamma activity at 60–90 Hz facilitates movement by seeking an inverse correlation between gamma activity and motor impairment.
It has been suggested that any prokinetic actions of gamma activity might relate to changes in arousal or motor vigor rather than to improvements in explicit motor processing (Kempf et al., 2009). This fits well with the view from studies in nonhuman primates suggesting that interactions between the basal ganglia and cortex are primarily modulatory in nature (Gatev and Wichmann, 2009). Our second aim was therefore to test the hypothesis that gamma activity is modulatory by comparing its reactivity during finger movements performed synchronously or asynchronously. Gamma increases of similar amplitude and duration at movement outset might be expected with modulatory processes such as arousal or motor vigor (Mazzoni et al., 2007), whereas more sustained increases in gamma activity might be predicted during sequential movements, if gamma is directly related to motor processing. In testing these hypotheses, we exploit recent methodological advances that permit the analysis of simultaneously recorded LFPs from deep neurosurgical targets and magnetoencephalographic (MEG) activity from the cortical mantle (Litvak et al., 2010; Hirschmann et al., 2011). The present report complements recent studies of cortico-subthalamic connectivity at rest in this (Litvak et al., 2011a) and other (Hirschmann et al., 2011) cohorts.
Materials and Methods
Patients and surgery.
Seventeen patients (55 ± 7 years old, six female, two left-handed) participated in the study. In one patient, the surgery was performed only on the right side. The study was approved by the joint ethics committee of the National Hospital for Neurology and Neurosurgery and the UCL Institute of Neurology, and the patients gave written informed consent before the study onset. Clinical details are given in Table 1. All patients were diagnosed with PD according the Queen Square Brain Bank criteria (Gibb and Lees, 1988). The indications, operative procedure, and beneficial clinical effects of STN stimulation have been described previously (Foltynie et al., 2011). Before surgery, the motor impairments of all patients were evaluated using part III of the Unified Parkinson's Disease Rating Scale (UPDRS) after omitting all dopaminergic medication overnight and after administration of at least 200 mg of levodopa. The patients were operated on after overnight withdrawal of levodopa medication and after dopamine agonists had been reduced and terminated during the 2 weeks before surgery.
Clinical features of patient cohort
The deep brain stimulation (DBS) electrodes were model 3389 (Medtronic) with four platinum–iridium cylindrical surfaces of 1.27 mm diameter, 1.5 mm length, and 2 mm center-to-center separation. The contacts were numbered 0 (lowermost, targeted to 2 mm below the center of the STN) to 3 (uppermost).
Surgical targeting of the DBS electrode was based on stereotactic magnetic resonance images (MRIs). Fast acquisition T2-weighted 2-mm-thick contiguous axial slices were acquired with a stereotactic Leksell Frame (Elekta). The STN [especially its medial border (Hariz et al., 2003)] was examined on the axial image containing the largest diameter of the ipsilateral red nucleus. The center of the STN was identified in a plane 0–1 mm behind the anterior border of the ipsilateral red nucleus (Bejjani et al., 2000). Cartesian coordinates of the target point were calculated using Framelink software (Medtronic). A double oblique trajectory was planned on reconstructed three-dimensional images to avoid sulci and ventricles. This surgical procedure has been described previously (Zrinzo et al., 2009; Foltynie et al., 2011). After implantation, electrodes were connected to an accessory kit, typically with both connectors being tunneled to the left temporoparietal area and externalized through the frontal region. No microelectrode recordings were made.
The locations of the electrodes were confirmed with immediate postoperative stereotactic imaging. Fast spin-echo T2-weighted 2-mm-thick contiguous axial slices were acquired with the Leksell frame still in situ. One patient was unable to tolerate a postoperative MRI and underwent stereotactic computed tomography scanning instead.
The patients were studied in the interval between DBS electrode implantation and subsequent connection to a subcutaneous stimulator between 2 and 7 d postoperatively.
Between 6 and 12 months postoperatively, the UPDRS assessment was repeated during DBS stimulation after overnight withdrawal of dopaminergic medication (Table 1, UPDRS OFF medication/ON DBS).
Experimental paradigm.
The experiment was divided into blocks of several minutes of recording. The patients could rest or fidget between the blocks. Each block comprised either rest or a movement task. During rest blocks, the subjects were instructed to remain still with their eyes open for 3 min (see Litvak et al., 2011a). During movement blocks, the subjects performed either simultaneous button presses with index, middle, and ring finger (SYN) or sequential button presses with index, ring, and then middle finger (SEQ), with either left or right hand (one kind of movement with the same hand within a block). The movements were self-paced. The subjects were instructed to move when they wanted but not to do it too frequently and to take ∼15 s between movements without counting silently. Feedback to the subjects was presented visually, using MATLAB (MathWorks) and a custom script based on Cogent (http://www.vislab.ucl.ac.uk/cogent.php). This script monitored the movement times and displayed messages on the screen when the intermovement interval was shorter than 12 s or a movement sequence was incomplete. When performing correctly, the subjects did not get any feedback and the screen just showed a fixation cross. In the case of incorrect performance, the script waited to collect another movement so that it ran until eight correctly performed movements were collected or at most for a total of 7 min. The subjects could usually complete a movement block in 3–4 min. A neurologist was present in the magnetically shielded room during the experiments to monitor the patients and performance of the task. The order of the conditions was pseudorandomized separately for blocks 1–9 and blocks 10–18, so that each half would contain one rest block and two movement blocks of each type. In all experiments, at least nine blocks were recorded, but only in two experiments did the patients successfully perform all 18 blocks. A single recording session lasted ∼1.5–2.5 h (including preparation).
The whole experiment was repeated twice: after overnight withdrawal of dopaminergic medication (OFF drug) and after the patients took ≥200 mg of levodopa (ON drug). The order of these drug conditions was counterbalanced over patients. Twelve patients were able to complete both experiments, three were only recorded ON drug, and two patients could not tolerate levodopa and were therefore only recorded OFF drug. Note that overnight withdrawal of dopaminergic medication is unlikely to lead to complete absence of dopaminergic input to the basal ganglia, because some dopaminergic effects may persist and some medications, particularly dopamine agonists, are long acting. Hence, the difference between the OFF and ON states is relative rather than absolute.
The two movement patterns differed in several regards, although the same fingers were used to enact the same key presses. In the SYN task, all three fingers simultaneously depressed three keys, whereas in the SEQ task, only one finger at a time was used to depress one key. Consequently, the conditions differed in the initial parameterization of force, which was greater in the SYN task. Nevertheless, the same three keys were eventually depressed in the SEQ task, which therefore went on for longer than synchronous button presses. Accordingly, activities directly related to motor processing might be expected to differ between tasks, reflecting their differential force profiles and durations. Moreover, the sequential movement was chosen to elicit motor-related processing in the basal ganglia, which are thought to be engaged by sequential movements (Benecke et al., 1987).
LFP–MEG recordings.
MEG recordings were obtained with a 275-channel system (CTF/VSM MedTech). Simultaneously, STN-LFP, electro-oculographic (EOG), and electromyographic (EMG) signals were recorded using an integrated EEG system (1× amplification, ±125 mV dynamic range) and high-pass filtered (in hardware) >1 Hz to avoid saturation of the amplifiers attributable to DC offsets. All signals were low-pass filtered in hardware <600 Hz. The data were sampled at 2400 Hz and stored to disk. Four intracranial LFP channels were recorded from each electrode, on each side, and referenced to a cephalic reference (forehead for the first two patients, right mastoid for the rest). LFP recordings were converted offline to a bipolar montage between adjacent contacts (three bipolar channels per hemisphere; 01, 12, and 23) to limit the effects of volume conduction from distant sources. EMG data were recorded from right and left first dorsal interosseous muscles with a reference at the muscle tendon. Button presses were also recorded in all subjects.
Head location was monitored using three head-position indicator (HPI) coils attached to the subject's nasion and pre-auricular points. For all but the first subject, we used continuous head localization and recorded the head locations throughout the experiment. Loss of head tracking occurred intermittently in some patients, possibly as a result of metal artifacts disrupting the head-tracking function of the MEG sensors. During offline processing, we compared the instantaneous distances between HPI coils with the distances based on the robust average (Holland and Welsch, 1977) of locations across the whole continuous recording. Time frames in which discrepancies were detected were discarded and replaced with linear interpolation based on the other time frames. This method works well when the tracking is valid for more than half of the recording, which was the case for all recordings reported here.
Data analyses.
In our previous study (Litvak et al., 2010), we detailed the challenges posed by the presence of metal artifacts in the MEG and our use of beam forming to suppress these artifacts. Here, we recap the main points and then focus on the details of the analysis specific to the present study.
Briefly, we used beamforming methods to localize the hand area of the primary motor cortex (M1) contralateral to the movement individually in each patient and extract virtual electrode time series from these locations. We then computed time–frequency images of M1 and STN power and M1–STN coherence around the times of the button presses by the contralateral hand and subjected these images to statistical analysis (see below, Statistical analysis). In the following sections, we describe the technical details of the beamforming and spectral analyses.
The data were analyzed using custom MATLAB scripts based on SPM8 (Litvak et al., 2011b) and Fieldtrip (Oostenveld et al., 2011) toolboxes (the Fieldtrip code we used is included in the SPM8 distribution).
Cortical source reconstruction.
For the purposes of M1 localization, we used a multivariate beamforming approach (Soto et al., 2009; Barnes et al., 2011). The motivation for using this approach was the ability to simultaneously examine a number of covarying spectral changes over a wide frequency range (5–100 Hz). When the induced responses at high (gamma) and low (alpha and beta) frequencies are localized separately (Jurkiewicz et al., 2006; Cheyne et al., 2008; Litvak et al., 2010), one is likely to get two separate peaks, and it is not clear which of them to use for virtual electrode placement. If gamma and alpha are covariant, then the optimal test is a multivariate one, which produces a statistic of the amount of variance (in the design) explained by a weighted linear combination of these two (or any N) spectral features. In this case, we used a feature space corresponding to the first five principal components of the power spectra at each virtual electrode location. For each patient and hemisphere, all contralateral movement trials were used together. The sensor locations were recomputed based on continuous head localization data (if available) to correspond to the mean location across the trials included in the analysis. Source time series were estimated with a linearly constrained minimum variance (LCMV) beamformer (Van Veen et al., 1997). The multivariate analysis contrasted two time windows: 0 to 1 s (activation) and −6 to −5 s (baseline) relative to the button press. The activation window was chosen as the period containing strong movement-related changes in power in both gamma and beta bands (see Results). The baseline window—matching the length of the activation window—was placed at the end of the baseline period used for time–frequency analysis (see below). χ2 statistical images were computed on a grid defined in Montreal Neurological Institute (MNI) space with spacing of 5 mm and restricted to the points within the inner skull boundary. Values on the grid were then interpolated using linear interpolation to produce volumetric images with 2 mm resolution. The resulting images were further smoothed with a 5 mm isotropic Gaussian kernel to remove spurious local maxima. We chose as individual M1 location the local maximum closest to the location defined in the literature (Mayka et al., 2006) [MNI coordinates: left M1, (−37, −18, 53); right M1, (37, −18, 53)].
Estimates of cortical source activity in these beamformer-identified M1 hand areas were examined in more detail using epochs from 10 s before to 10 s after button presses. Again, virtual electrode time series were estimated with LCMV beamformer (Van Veen et al., 1997). The source orientation was set in the anteroposterior direction to be approximately normal to the cortical surface in the central and precentral sulcus.
When computing covariances for beamforming, it is necessary to account for head movement. Ideally, one should use as much data as possible (Brookes et al., 2008). For instance, when looking at movement-related activity, it would be desirable to use all movement trials. However, when the head location varies across trials, the sources seen by a particular MEG sensor may differ and pooling over trials may then degrade covariance estimates. In experiments with healthy volunteers, this is a minor issue because head movements do not usually exceed 1 cm (which is of the order of the sensor coil radius in our MEG system). In our case, however, some patients found it difficult to keep still because of involuntary dyskinesias, especially when on dopaminergic medication. Thus, displacements of ∼3 cm between trials were not uncommon. We, therefore, computed beamformer filters for each trial separately. This is also the most generic approach accommodating varying numbers of trials across subjects. In all cases, the sensor positions were recomputed to correspond to the average head location for each trial. We also tested two alternative approaches—computing filters based on all the trials in an experiment and on all the trials in each block—and obtained qualitatively similar results. Furthermore, we tested a range of regularization values that constrain the beamformer estimates. This indicated that the minimum regularization value yielding consistent results was 0.01% of the signal variance (averaged over channels).
Preprocessing of virtual electrode and LFP data.
The virtual electrode channels derived from the MEG sometimes contained discontinuous jumps, whose origin could be traced to occasional resets of the sensor circuitry. These jumps were detected by thresholding the differences between adjacent samples. When a jump was detected, the values from 20 samples before to 20 samples after the jump were replaced by the median difference over this segment and the modified difference time series were summed again to produce the original time series with the jump corrected.
The corrected data were digitally filtered (1 Hz fifth-order high pass and fourth-order notch filters for 50 Hz and all harmonics up to 550 Hz, zero-phase Butterworth in all cases). Finally, before spectral analysis, the channel data were standardized by subtracting the mean and dividing by the standard deviation (SD) for each channel and trial separately. This ensured that all trials contributed equally to estimates of source activity.
Artifact suppression.
High-frequency modulations in power, particularly wide-band modulations, could result from artifacts in the forms of short “spikes” or “jumps” in the data. In our case, LFP, EMG, and EOG data were contaminated by occasional brief electrical discharges attributable to a grounding problem that could not be completely resolved at the hardware level.
To rule out high-frequency artifacts as an explanation for our results, the pre-processed LFP and virtual electrode data and their first derivative (difference) time series were examined for the presence of outliers by thresholding. The thresholds were set such that they could separate the artifacts from the remaining data. These thresholds were 5 SDs for the original data, 4 SDs for STN–LFP derivative, and 1.5 SDs for M1 derivative (Fig. 1A,C).
Artifact rejection criteria and validation of robust averaging. A, Fifty trials from one subject showing the spike artifacts in STN–LFP. The red lines indicate the rejection threshold that was used for STN–LFP data (5 SDs). B, Distribution of STN–LFP artifacts in trial time. The rate of artifact occurrence increased around the button press, peaking at ∼0.5%. If an artifact occurred within the boundaries marked by the red lines, the corresponding trial was excluded from analysis of power. C, Difference time series for 50 trials of M1 virtual electrode data. The red lines indicate the rejection threshold that was used for these data (1.5 SDs). D, Distribution of M1 artifacts in trial time. The rate of artifact occurrence increased around the button press, peaking below 1%. If an artifact occurred within the boundaries marked by the red lines, the corresponding trial was excluded from analysis of power. E, Results of simulations aimed at validating the robust averaging method. The top row shows the results of time–frequency analysis for synthetic data contaminated with artifacts taken from the real data. The second row shows similar analysis on simulated data without artifacts, and the third row shows the results of analysis of contaminated data when using robust averaging. The simulated data were based on real data from all the subjects used in the actual data analysis and one experimental condition (synchronous button presses with the right-hand ON drug). To assess reproducibility of these results, the simulation was repeated with eight different sets of artifacts taken from all the eight experimental conditions. The results were summarized by averaging the power and coherence in 15–35 Hz range. The bottom row shows the mean and SD of these eight repetitions for contaminated signals (black), clean signals (red), and robust averaging results (green). The units in all plots are change from baseline in percentage.
Figure 1A shows an example of LFP recording containing spikes from a single experiment, and Figure 1C shows an example of M1 signal derivative also evidencing artifacts in some trials. Examining the timing of the artifacts relative to the button press showed that their occurrence rate was fairly uniform, with the exception of a short segment around the button press, in which the incidence of spikes increased (Fig. 1B,D). To preclude contribution of spikes to the button-press response, trials with artifacts in this period (−0.2 to 1.5 s) were excluded from the analysis of event-related power. The effects of the remaining spikes were suppressed by robust averaging (see below). Trial exclusion was performed separately for M1 and STN channels to minimize data loss. For the analysis of event-related coherence, we used all the trials. This was motivated by simulations showing that high-frequency artifacts in only one of the channels had a minor effect on coherence estimates and were further suppressed by our robust procedure for coherence computation (see below, Robust averaging).
The numbers of trials included in the analysis for each subject, hemisphere, and condition are reported in Table 2.
Details of hemispheres included in the analysis
Spectral analyses.
For efficient spectral estimation from a relatively small number of trials, we used multitaper spectral analysis (Thomson, 1982). This method is based on pre-multiplying the data with a series of tapers optimized for producing uncorrelated estimates of the spectrum in a given frequency band. This sacrifices some of the frequency resolution, in a controlled manner, to increase the signal-to-noise ratio. It does this by effectively multiplying the number of trials by the number of tapers used. We estimated the spectra between −8 and 8 s relative to the first button press of each trial, in overlapping windows of 400 ms (shifted by 50 ms). The frequency resolution was set to the inverse of the time window (2.5 Hz) for up to 25 Hz, then 0.1 times the frequency for 25 to 50 Hz, and then to a constant 5 Hz resolution. These settings resulted in a single taper being used for 2.5–30 Hz, two tapers for 32.5–42.5 Hz, and three tapers for ≥45 Hz. The resulting time–frequency images had no discontinuities thanks to the continuous frequency resolution function.
The time–frequency images were then averaged using robust averaging (Holland and Welsch, 1977; Wager et al., 2005; for details, see below), and percentage change time–frequency responses were obtained by normalizing to the baseline (−8 to −5 s) before button press.
Coherence was estimated using a similar spectral estimation procedure except that robust averaging was used during coherence computation and was applied separately to the absolute values of the cross-spectral density and to the power of the two sources. In the case of cross-spectra, the weights computed from the absolute values were then applied to the complex cross-spectra when computing the mean coherency over trials. Percentage changes in cross-spectral responses were computed as above.
For the purposes of data reduction, we only selected one of the three STN bipolar channels on each side. The initial selection was done based on the greatest movement-related beta event-related desynchronization (ERD)/event-related synchronization (ERS). We then also examined the average baseline-corrected spectra around the button press to ensure that no other spectral features were completely omitted by selecting this particular channel and changed the selection to a more representative channel if necessary. The same channel was also used for coherence computation. This approach is based on the fact that the highest beta activity in the subthalamic region localizes to the STN, particularly its dorsal “motor” part, as judged by microelectrode recordings of neuronal discharges and LFPs in parkinsonian models (Sharott et al., 2005; Mallet et al., 2008; Degos et al., 2009) and patients (Kühn et al., 2005; Chen et al., 2006; Trottenberg et al., 2006; Weinberger et al., 2006; Moran et al., 2008; Yoshida et al., 2010; Zaidel et al., 2010).
Robust averaging.
Robust averaging is a special and simple case of the robust general linear model (Holland and Welsch, 1977; Wager et al., 2005). The idea is that, for each time–frequency bin, the distribution of values over trials is considered and the outliers are down-weighted when computing the average. This makes it possible to suppress artifacts restricted to narrow time and frequency ranges without rejecting whole trials. Moreover, a clean average can be computed with no clean trials, provided that artifacts do not overlap systematically. This was the case for our data, except in some trials around the button press, in which we excluded artifacts by thresholding (see above).
To verify the validity of robust averaging for our data, we applied it to simulated data combined with spike artifacts from the STN–LFP recordings. The simulated data were derived from the original data; STN–LFP data for each trial were replaced by a simulated signal generated as follows. The amplitude of bandpass-filtered (15–35 Hz) noise was modulated to create an increase in amplitude by 50% around time 0. The amplitude modulation was Gaussian, with SD of 250 ms. This narrow-band noise was added to white noise (in equal proportions). The resulting signal was standardized by subtracting the mean and dividing by the SD for each channel and trial separately to create the “clean” reference LFP signal (henceforth “clean signal”). Contaminated LFP signal (henceforth “contaminated signal”) was created by adding the spikes from the original pre-processed LFP data from the corresponding trial to the clean signal. The spikes were separated from physiological signal by setting the values with amplitude <5 SDs to 0. For the purposes of estimating coherence, the clean signal was added to real M1 virtual electrode signal to create a signal coherent with the clean and the contaminated signals (henceforth “coherent signal”). We then subjected the synthetic signals to the same time–frequency analysis as real data with and without the use of robust averaging for power and coherence computation. The simulated data were generated using real data from all subjects and one experimental condition (synchronous button presses with the right hand ON drug). To assess reproducibility of the results, the simulation was repeated with eight different sets of artifacts taken from all eight experimental conditions. For the purposes of the simulation, we used all the original trials, including the trials with artifacts around the button press, which makes the results for power conservative with respect to the actual analysis.
Figure 1E shows the simulation results. For power estimates, spike artifacts strongly affect the results as would be expected. However, robust averaging enables one to recover the simulated power increase. For coherence estimates, the simulated effect can still be seen even when contaminated data are used. The reason for this is that the power associated with a spike in one channel only affects the relative amount of shared power at any given frequency. However, robust averaging is also advantageous in this case, because it reduces the variance in coherence estimates (Fig. 1E, bottom row), thereby increasing the sensitivity of the statistical analysis.
Based on the results of these simulations, trials with artifacts in the window from −0.2 to 1.5 s around the button press were excluded from the analysis of power but retained for the analysis of coherence (except in Figs. 7 and 9).
Excluding data with high frequencies in the evoked response.
To further ensure that the phenomena we report pertain to induced rather than evoked responses, we computed the averages of both the movement-related LFP and virtual electrode data in the time domain and performed time–frequency analysis of these evoked responses using the same settings as for single trials. Even after exclusion of trials containing artifacts, for some hemispheres, we found high-frequency activity around the button press. All of these cases were from M1 virtual electrode data and not from STN–LFPs. Not all of this activity was clearly artifactual, but because this was a likely explanation, we excluded all hemispheres in which such activity was found from our analysis of M1 power. Eight hemispheres were excluded, and, for one additional subject (two hemispheres), only the data recorded ON medication were excluded. Table 2 details the data included in the analysis.
Statistical analyses.
To make inferences about the spectral responses and effects of task and drug, we used statistical parametric mapping (as implemented in SPM8). This treats the time–frequency matrices as images and allowed the identification of regions of time–frequency space showing significant effects over subjects, while controlling for the implicit multiple comparisons using random field theory (Kilner et al., 2005; Kilner and Friston, 2010).
The results of time–frequency analysis were exported to Neuroimaging Informatics Technology Initiative format and smoothed with a Gaussian smoothing kernel with full-width half-maximum of 7.5 Hz by 500 ms for frequencies 0–100 Hz and 25 Hz by 500 ms for frequencies 100–600 Hz.
All the reported findings are significant with familywise error (FWE) correction at the cluster level (p < 0.01 corrected, cluster forming threshold p < 0.05 uncorrected).
Significant features of the mean responses were determined by subjecting mean images across conditions for power and coherence to a single-sample t test across subjects. t tests were performed separately for dominant and nondominant hemispheres so that data from each subject were only used once in each test.
To test for the effects of experimental conditions, we performed repeated-measures ANOVA, treating the experiment as a 2 × 2 factorial design with the factors “task” (SYN vs SEQ) and “drug” (ON vs OFF). Here, data from dominant and nondominant hemispheres were pooled and modeled as dependent to account for the resulting correlations in the error.
To test for correlation with clinical scores, we added to the ANOVA described above a regressor with contralateral hemibody bradykinesia-rigidity scores. This analysis was done separately because the effect of drug is highly correlated with these scores. The scores comprised the sum of items 22–26 of the UPDRS part III score. Motor examination was performed preoperatively, after withdrawal from medication overnight, in a practically defined off state (so that patients had their last antiparkinsonian medication 9–12 h before testing). Motor examination was repeated on the same day 1 h after their usual antiparkinsonian treatment, provided that the levodopa dose was ≥200 mg. When this was not the case, the patient's standard antiparkinsonian medication was replaced by a single dose of 200 mg of levodopa. The contralateral hemibody bradykinesia-rigidity score was selected as the clinical regressor because, unlike contralateral tremor score, this has been found repeatedly to correlate with LFP (beta-band) activity in the subthalamic region (Kühn et al., 2006, 2009; Ray et al., 2008; Chen et al., 2010).
Localization of gamma coherence.
For localization of gamma coherence in several example subjects, source localization was performed using the dynamic imaging of coherent sources (DICS) beamforming method (Gross et al., 2001; Litvak et al., 2010). Before beamforming, trials with artifacts in the LFP data were excluded by thresholding the (nonstandardized) LFP at 20 μV. Beamformer settings were as for the multivariate beamformer (see above).
Directionality analyses.
To determine the directionality of functional coupling, we compared coherence with nonparametric variant of Granger causality (Brovelli et al., 2004). Both measures were computed based on multitaper spectral analysis of two time windows: baseline (−6 to −5 s) and movement (−0.5 to +0.5 s) relative to the button press with spectral resolution of 5 Hz. The directionality analysis did not use robust averaging, and trials in which deflections exceeding 5 SDs occurred in the analysis window were excluded from the analysis.
Results
Localization of M1
Seventeen patients participated in the study. In one of them, LFP recording was only possible from one hemisphere. Thus, localization of M1 using the multivariate beamformer was performed in 33 hemispheres. Locations of individual beamformer peaks for each hemisphere are reported in Table 2. The distance from individual beamformer peaks transformed to MNI space to the location of hand area of M1 was 18.6 ± 6.1 mm (mean ± SD; 6.8 mm minimum, 33 mm maximum).
Responses in M1 and STN induced by voluntary movement
Figure 2 shows the power changes induced by button presses with the contralateral hand (relative to baseline; mean over all conditions). The top row of Figure 2A shows (unthresholded) time–frequency images, and the bottom row reports the corresponding significant clusters.
Power responses induced by the button press in contralateral M1 and STN. A, Average time–frequency images. Induced responses were baseline corrected (baseline −8 to −5 s) and averaged across conditions. Separate t tests were performed for dominant and nondominant hemispheres and for 0–100 and 100–600 Hz ranges (each condition has 2 spectrograms). The top row shows unthresholded mean time–frequency images relating to M1 (leftmost 4 panels) and STN power change (rightmost 4 panels), and the bottom row reports the corresponding significant positive (white) and negative (black) clusters (p < 0.01, cluster-level FWE correction). B, Induced responses of individual hemispheres. The individual responses whose averages are presented in A were averaged between 0 and 1 s relative to the button press. Dominant and nondominant hemispheres are presented together. Note the presence of clear gamma activity in the 40–90 Hz range in many of the individual hemispheres as well as 300–400 Hz activity clearly present in four hemispheres (for details see Table 2).
Separate t tests were performed for 0–100 and 100–600 Hz ranges. A significant ERD in the beta band and an ERS in a broad (40–250 Hz for M1, 40–600 Hz for STN) gamma band occurred at the time of the movement onset in both the M1 and STN contralateral to the movement. In both dominant and nondominant STNs, there was a distinct peak in the ERS between 300 and 400 Hz.
Figure 2B shows induced responses for individual hemispheres averaged between 0 and 1 s relative to the button press. In M1, most hemispheres displayed a distinct peak in the 50–100 Hz range.
In the STN, discrete peaks in individual subjects were observed under 100 Hz and in the 300–400 Hz range. The 300–400 Hz activity in response to movement was observed in four hemispheres (three subjects) detailed in Table 2. In both M1 and STN, patients with and without a given spectral profile did not differ in terms of disease duration, disease severity, or dominance of tremor or akinesia rigidity (data not shown).
Effects of experimental condition on induced responses
The effect of drug at M1 was manifest as a wide-band increase in induced activity around first button press (Fig. 3). The response in STN was also wide band but more extended in time and with distinct peaks around 10, 70, and 300 Hz.
Effect of drug on power responses induced by button press in the contralateral M1 and STN. The top row shows the unthresholded contrast images corresponding to the effect of drug in the ANOVA, and the bottom row reports the corresponding significant clusters (p < 0.01, cluster-level FWE correction).
The effect of task on M1 gamma was observed in two different time windows with opposite directions (Fig. 4A). Figure 4B shows the induced responses in the 0–100 Hz range for each task separately (averaged over drug conditions) to aid interpretation of the difference image. Around the time of the first button press, the broad gamma band ERS was reduced during sequential relative to synchronous button presses, but this effect did not reach significance above 50 Hz. The greater gamma activity time-locked to the initial button press in the SYN task may have related directly or indirectly to the greater force required at this point in time in this task. This potentially force-related processing was seen in M1 and the STN but was more marked in M1. Note, however, that gamma power increased relative to the baseline in both cases, and the negative effect emerges because we are considering the difference image between two conditions. After the initial button press, the gamma response in the sequential task was significantly stronger than with synchronous button press. The beta ERS during sequential button press was reduced relative to that during synchronous button presses. In the STN, similar differences were observed but they were much weaker and none reached significance. Notably, the gamma response in the STN is prolonged in both SYN and SEQ tasks in contrast to the brief gamma response in M1 for the SYN task. This difference cannot be explained by exclusion of some of the hemispheres from the M1 power analysis, because excluding the same hemispheres from the STN analysis had no effect on the results (data not shown). There were no significant interactions between drug and task effects.
Effect of task on power responses induced by button press in the contralateral M1 and STN. A, The top row shows the unthresholded contrast images [in pairs of low (0–100 Hz) and high (100–600 Hz) frequency spectrograms] corresponding to the effect of task in the ANOVA, and the bottom row reports the corresponding significant clusters (p < 0.01, cluster-level FWE correction). B, Averaged power responses (in the 0–100 Hz band only) in M1 and STN for each of the two tasks separately.
Reactivity patterns of M1–STN coherence differ from those of local synchronization
Figure 5 shows baseline corrected M1–STN coherence contralateral to movement. There was an increase in coherence over a narrow gamma band from 65 to 80 Hz at the time of the first button press and shortly thereafter. This effect only reached significance for the dominant hemispheres. For both dominant and nondominant hemispheres, there was also a significant coherence increase in the beta band starting at ∼2 s after the button press, which might be related to the beta band power increase at this time.
Contralateral M1–STN coherence responses induced by the button press. Coherence images were baseline corrected (baseline −8 to −5 s) and averaged across conditions. Separate t tests were performed for dominant and nondominant hemispheres and for 0–100 and 100–600 Hz ranges. Only 0–100 Hz results are shown because for 100–600 Hz there were no significant effects. The top row shows unthresholded mean time–frequency images, and the bottom row reports the corresponding significant clusters (p < 0.01, cluster-level FWE correction).
Figure 6 shows the effects of experimental condition on coherence. The increase in narrow gamma (65–80 Hz) band coherence at the time of button press was further increased by treatment with dopaminergic medication. For five hemispheres (three subjects), this effect could be seen as distinct peak in the difference between ON and OFF drug spectra (Fig. 6B; for details, see Table 2). To confirm that this effect was specific to M1, we localized the gamma coherence for these hemispheres with DICS. For the purposes of localization, we used combined data from both tasks ON drug, time window of −0.5 to 0.5 s relative to the button press and frequency windows centered around the individual coherence peaks. In three of these hemispheres, the topography of the narrow gamma band coherence with depth LFPs was centered on ipsilateral M1, but in two (from the same subject), coherence was centered a few centimeters anterolateral to M1. The peak values are given in Table 2 (footnotes), and an example is shown in Figure 6C. There were no significant effects of task and no significant interaction of task and drug.
Effects of experimental condition on M1–STN coherence responses. A, The top row shows the unthresholded contrast of M1–STN coherence images corresponding to the effects of drug and task in the ANOVA, and the bottom row reports the corresponding significant clusters (p < 0.01, cluster-level FWE correction). B, Effect of drug on coherence responses in individual hemispheres. Differences between ON and OFF drug coherence images (collapsed across the task factor) were averaged between −0.5 and 0.5 s relative to the button press. The red line corresponds to the localization results in C. C, Localization of gamma coherence in an individual subject (subject 11 left hemisphere). DICS beamformer was applied to button presses with the contralateral hand ON drug pooled across tasks. The time range for the analysis was −0.5 to 0.5 s relative to the button press, and the frequency range was 65–85 Hz. The image was transformed to MNI template space and overlaid on the template structural image. The peak coherence was observed at MNI coordinates (−30, −4, 64).
To further verify that narrow band gamma coherence in our data was not caused by spike artifacts, we repeated the coherence computation for the data ON drug from the five hemispheres showing clear coherence peaks after excluding trials with artifacts in the window −1 to 1 s relative to the button press in either M1 or STN signal. The criteria for exclusion were the same as for power analysis (see Materials and Methods and Fig. 1). The resulting raw data and coherence spectra (without baseline correction) are shown in Figure 7. The coherence peaks in the gamma range were preserved and therefore could not be explained by artifacts. Note that we could not apply the same exclusion criteria to all the data because this would result in too few trials for some subject–condition combinations.
Narrow band coherence peaks in the gamma range were not caused by artifacts. Coherence computation was repeated for the data ON drug from the five hemispheres showing clear gamma coherence peaks after excluding trials with artifacts in the window −1 to 1 s relative to the button press in either M1 or STN signal. The criteria for exclusion were the same as for power analysis (see Materials and Methods and Fig. 1). Time–frequency decomposition was the same as for Figures 5 and 6, and the coherence values were averaged between −0.5 and 0.5 s relative to the button press. Raw data for all the trials are shown as well as the coherence spectra. The subject identification numbers and the trial numbers are detailed on the right.
Correlations with motor state
To assess the possible clinical relevance of the spectral changes described above, we performed ANCOVA with the contralateral hemibody bradykinesia-rigidity scores as an additional regressor (Fig. 8). The sign of the scores was flipped so that positive correlations mean that drug-induced increase in reactivity of power or coherence during movement correlates with drug-induced behavioral improvement. No significant correlations were found with M1 power. Significantly correlated clusters of STN power and M1–STN coherence were identified in the 60–90 Hz band around the time of, and just after, the initial button press. These were the same times and frequencies for which there were also significant effects of drug. There was also a significant positive correlation with the 300 Hz power in the STN after the button press. Finally, there was a negative correlation with M1–STN coherence between 15 and 45 Hz at the time of the button press. The same features were observed with pure bradykinesia hemibody scores as a regressor (data not shown). To demonstrate that the correlations with STN power and M1–STN coherence in the 60–90 Hz range were not driven by outliers, we also plotted the relations between these effects and clinical improvement for individual hemispheres (Fig. 8B).
Correlations between drug-induced changes in power and coherence responses with drug-induced changes in clinical scores. Repeated measures ANCOVA was performed with contralateral hemibody bradykinesia-rigidity scores as an additional regressor. The sign of the scores was changed to negative so that positive correlations would correspond to clinical improvement. A, The top row shows unthresholded contrast images corresponding to the effect of clinical score in the ANCOVA (the units are percentage change from baseline per unit change in motor UPDRS), and the bottom row reports the corresponding significant positive (white) and negative (black) clusters (p < 0.01, cluster-level FWE correction). B, Relationship between the effects of drug in STN power and coherence responses in the 60–90 Hz range and clinical improvement. Power and coherence responses were averaged over −0.5 to 0.5 s window. The solid lines show linear fit to the data (for STN power, r2 = 0.27, p = 0.01; for coherence, r2 = 0.42, p = 0.001).
Directionality analyses
To determine the directionality of the interactions between M1 and STN, we compared coherence computed for individual hemispheres with nonparametric Granger causality analysis. This analysis was done for two 1 s windows: baseline (−6 to −5 s) and movement (−0.5 to +0.5 s) relative to the button press. The results are shown in Figure 9. The baseline coherence spectrum is dominated by beta coherence, with Granger causality indicating a predominant M1 to STN direction as confirmed by a paired t test between averages of Granger causality for the two directions in the 15–30 Hz range (p < 0.02). The movement window is characterized by the presence of 60–90 Hz peaks with STN to M1 directionality. Paired t test between averages of Granger causality for the two directions in the 60–90 Hz range showed this effect only as a trend (p < 0.1) because of the small number of hemispheres with clear directional peaks in the gamma range.
Directionality of M1–STN coupling. Coherence and nonparametric Granger causality were computed for individual hemispheres (the data were combined across conditions). This analysis was done for two 1 s windows: baseline, −6 to −5 s, and movement, −0.5 to +0.5 s relative to the button press. Note that, in the baseline period, coherence is present in the beta band with the predominant direction being from M1 to STN, whereas in the movement period, coherence is present in the gamma band with the predominant direction being from STN to M1. In both cases, there is clear correspondence between coherence and Granger causality results.
Discussion
In this study of patients with PD, we focused on local and long-range synchronization in the gamma frequency band. Such activity was polymorphic with several features, distinguished by their extent and reactivity. These consisted of discrete peaks in STN power at 60–90 and 300–400 Hz and a wide-band spectral shift extending to 600 Hz. All increased with movement and levodopa treatment. Only STN gamma activity at 60–90 Hz was coherent with activity in M1, which it tended to drive. We obtained some evidence that this coherence could be specific to the motor areas ipsilateral to the STN. However, because our paradigm was not optimal for collecting sufficient data to localize coherence, additional studies are necessary to confirm this finding. Drug effects on both local and distant synchronization at 60–90 Hz correlated with the degree of improvement in contralateral bradykinesia, as did local STN activity at 300–400 Hz. Although correlation does not necessarily imply causation, this finding does provide support for the hypothesis that such high-frequency activities are prokinetic in function. Thus, greater increases in power and coherence are associated with less clinical impairment, in contrast to the opposite findings in the beta frequency band made here and previously (Doyle et al., 2005). These relationships were significant despite the long interval between clinical assessments and recordings and the possible confound of postoperative “stun” or microlesional effects (Chen et al., 2006).
Previously, the only evidence in favor of a “prokinetic” action of the 60–90 Hz activity was its increase with voluntary movement and after treatment with levodopa (Brown et al., 2001; Cassidy et al., 2002; Alegre et al., 2005; Fogelson et al., 2005; Alonso-Frech et al., 2006; Devos et al., 2006; Androulidakis et al., 2007) and a trend toward a correlation between treatment-related increases in background 60–90 Hz power in the STN and clinical improvement (Kühn et al., 2006; Trottenberg et al., 2006). Our findings are all the more remarkable given that high-frequency synchronization has primarily been observed in short-range corticocortical interactions and usually in the context of perception (Fries, 2009). Although several previous studies reported long-range gamma synchronization (Williams et al., 2002; Schoffelen et al., 2005; Lalo et al., 2008), it is still often presumed that synchronization over long distances is unlikely to occur at high frequencies (Kopell et al., 2000; Ray and Maunsell, 2010). A previous study investigating the coupling between the STN–LFP and cerebral cortex during phasic movements in PD patients only sampled EEG from a limited number of sites, obviating even approximate cortical localization, and averaged movement periods with intervening periods without movement in a block design (Lalo et al., 2008).
Interestingly, the 60–90 Hz activity in STN and corresponding coherence, unlike activity at similar frequencies in M1, was unaffected by the type of movement performed. Thus, the increase in gamma activity was little different in terms of its size or duration between synchronous and sequential finger movements. The stereotyped nature of the induced gamma activity across tasks suggests that it might underlie some modulatory feature common to both motor acts, such as phasic arousal or attention (Kempf et al., 2009). The correlation with general motor state would be consistent with a modulatory function, as would the observation that the 60–90 Hz activity was not necessary for movement insofar as it was not always observed. Such 60–90 Hz activity has been identified in and between several sites in the basal ganglia–cortical loop across several disease states, suggesting that it may be primarily physiological rather than pathological (Brown et al., 2001; Cassidy et al., 2002; Kempf et al., 2009). It is modulated by movement and, importantly, varies over the sleep–wake cycle (Brown et al., 2001; Kempf et al., 2009). Moreover, it is enhanced by startle-eliciting stimuli (Kempf et al., 2009), consistent with modulation by arousal state.
The notion that basal ganglia input to motor cortical areas may be related to arousal is not new: Hassler (1978) first proposed that the brainstem reticular activating system impacted on the basal ganglia and thereby motor cortical areas, through the thalamus. Our studies support this hypothesis and suggest that 60–90 Hz activity may mediate this effect. Accordingly, we have demonstrated a phasic (direct or indirect) drive from the STN to the motor cortex in the gamma band around the time of movement. Consistent with this, Williams et al. (2002) showed that STN activity in this frequency band led coupled activity in the cerebral cortex in parkinsonian patients with prominent gamma activity at rest.
In line with previous studies, we also found a spectral peak at ∼300 Hz, which increased after levodopa and correlated with contralateral motor state (Foffani et al., 2003; Kane et al., 2009; López-Azcárate et al., 2010; Ozkurt et al., 2011). There was no evidence for STN–M1 coherence at these very high frequencies, suggesting that they are features of local processing rather than coupled activities in the distributed subcortico-cortical network.
The other gamma band feature was a movement-induced wide-band spectral power increase. This was evident in both M1 and STN. In the latter, it could occur with or without more discrete peaks at 60–90 and 300–400 Hz. The wide-band power change was promoted at both levels by treatment with levodopa. Such high-frequency changes in spectral power have been reported previously and hypothesized to be an LFP correlate of population firing rate (Manning et al., 2009; Miller et al., 2009; Ray and Maunsell, 2010).
The different forms of gamma activity were not evident in all our patients. This might reflect phenotypic variation or disease severity, although no clear association was apparent. Additional studies will be necessary to explore these potential sources of variation in a larger patient cohort. Alternatively, the variability in spectra between patients might be attributable to targeting variance, sampling error given the rather focal source of the gamma activity (Trottenberg et al., 2006), postoperative stun effects (Chen et al., 2006), or insufficient doses of levodopa in some cases. Note that all these activities were identified through their phasic relationship to movement, explaining why previous analyses of resting MEG data failed to identify clear cortical sources of gamma activity that were coherent with STN activity (Hirschmann et al., 2011; Litvak et al., 2011a).
The coherence between some STN and M1 activities raises the possibility of volume conduction. Many arguments have been put forward against this (Brown and Williams, 2005): most convincing is that the discharge of neurons in the STN tends to be locked to those oscillations in the STN–LFP that are coherent with cortex (Levy et al., 2002; Kühn et al., 2005; Trottenberg et al., 2006; Weinberger et al., 2006). It is also worth commenting that the patterns of oscillatory responses were not identical between M1 and STN. This is additional evidence against far-field volume conduction effects being picked up by our subthalamic electrode and points to either functional specialization of the two sites or their differential involvement in the disease process. Thus, task-related local processing was more prominent in M1 than STN, as evidenced by the long duration of broad-band gamma power increases in M1 but not in STN in sequential movements.
It should be noted that coherence is only a measure of linear coupling, and its absence does not rule out the possibility of nonlinear coupling, which can be much more diverse in its mechanisms and therefore more difficult to detect (Jensen and Colgin, 2007). Here we focused on linear coupling, and examining other possibilities will be part of our future work.
Another important issue is the extent to which the effects of treatment and movement could be attributable to changes in sensory re-afference secondary to changes in movement patterns (e.g., reduced bradykinesia). Although this may contribute to differences in responses after movement onset (Muthukumaraswamy, 2010), it cannot explain the altered changes in local power (Kühn et al., 2004; Loukas and Brown, 2004; Doyle et al., 2005; Williams et al., 2005; Androulidakis et al., 2007; Kempf et al., 2007; Alegre et al., 2010) and cortico-STN coherence (Cassidy et al., 2002), which start before the movement in both the beta and 60–90 Hz bands. Consistent with this, reactivity of STN–LFP power and cortico-STN coherence in the beta band still occurs in imagined movements (Kühn et al., 2006), during movement observation (Alegre et al., 2010), and is greatly attenuated during passive movements (Cassidy et al., 2002). Thus, at least some of the changes in spectral responses are likely to have been primary and not an epiphenomenon of re-afference, although the latter remains a possibility with respect to the movement-induced wide-band spectral shift and activity centered on 300 Hz.
In conclusion, we have demonstrated several movement-related increases in the gamma band. Of these, the 60–90 Hz band power in STN and corresponding coherence with M1 correlates with contralateral bradykinesia rigidity in patients with PD, supporting the prokinetic nature of this form of synchronization. That said, this synchronization was independent of the type of movement performed and may therefore underpin a general feature of motor control.
Footnotes
This study was funded by Marie Curie Intra European Fellowship MEIF-CT-2006-038858 and a Coleman-Cohen fellowship of the British Technion Society to V.L. A.J. is supported by a training fellowship from Parkinson's UK, and M.I.H., L.Z., T.F., and P.L. are supported by the UK Parkinson's Appeal. P.B. is supported by the Medical Research Council, Rosetrees Trust, and the National Institute of Health Research Biomedical Research Centre, Oxford. V.L., K.F., and A.E. are supported by the Wellcome Trust. We thank David Bradbury, Janice Glensman, and Zoe Chen for their assistance conducting the experiments.
- Correspondence should be addressed to Dr. Vladimir Litvak, Wellcome Trust Centre for Neuroimaging, 12 Queen Square, London, WC1N 3BG, UK. v.litvak{at}ucl.ac.uk
