Abstract
In Parkinson's disease (PD), striatal dopamine denervation results in a cascade of abnormalities in the single-unit activity of downstream basal ganglia nuclei that include increased firing rate, altered firing patterns, and increased oscillatory activity. However, the effects of these abnormalities on cortical function are poorly understood. Here, in humans undergoing deep brain stimulator implantation surgery, we use the novel technique of subdural electrocorticography in combination with subthalamic nucleus (STN) single-unit recording to study basal ganglia–cortex interactions at the millisecond time scale. We show that in patients with PD, STN spiking is synchronized with primary motor cortex (M1) local field potentials in two distinct patterns: first, STN spikes are phase-synchronized with M1 rhythms in the theta, alpha, or beta (4–30 Hz) bands. Second, STN spikes are synchronized with M1 gamma activity over a broad spectral range (50–200 Hz). The amplitude of STN spike-synchronized gamma activity in M1 is itself rhythmically modulated by the phase of a lower-frequency rhythm (phase-amplitude coupling), such that “waves” of phase-synchronized gamma activity precede the occurrence of STN spikes. We show the disease specificity of these phenomena in PD, by comparison with STN-M1 paired recordings performed in a group of patients with a different disorder, primary craniocervical dystonia. Our findings support a model of the basal ganglia-thalamocortical loop in PD in which gamma activity in primary motor cortex, modulated by the phase of low-frequency rhythms, drives STN unit discharge.
Introduction
Much work has been done to define the electrophysiology of the basal ganglia nuclei in humans with movement disorders. In Parkinson's disease (PD) patients off dopaminergic medications, subthalamic nucleus (STN) neurons have increased neuronal firing rates (Hutchison et al., 1998; Steigerwald et al., 2008; Schrock et al., 2009), oscillatory firing patterns in the theta (4–8 Hz), alpha (8–12 Hz), and beta (12–30 Hz) ranges (Rodriguez-Oroz et al., 2001; Levy et al., 2002b), and exaggerated synchronization to neighboring STN units and to STN local field potentials (LFPs) in the beta range (Levy et al., 2002b; Kuhn et al., 2005; Weinberger et al., 2006; Moran et al., 2008). Excessive basal ganglia neuronal spike synchronization is also well established in animal models of parkinsonism (Nini et al., 1995; Goldberg et al., 2004; Costa et al., 2006).
The role of the motor cortex in parkinsonian pathophysiology is less clear. Scalp electroencephalography and invasive LFP studies in humans without movement disorders show that beta frequency oscillations are normally present in motor cortex (Pfurtscheller et al., 1997; Crone et al., 1998; Miller et al., 2007). In rodent models, cortical beta rhythms are excessively synchronized to basal ganglia spike discharge (Mallet et al., 2008; Brazhnik et al., 2012), but it is not known whether this is true in humans with PD. Further, interactions between basal ganglia spiking and higher frequency cortical activity have not been explored. The broad range of gamma frequencies between 50 and 200 Hz, often referred to as “broadband gamma” (Ossandon et al., 2011), is probably a surrogate measure of underlying asynchronous spiking (Manning et al., 2009), and task-related changes in broadband gamma are thought to reflect localized cortical function (Lachaux et al., 2012). Further, broadband gamma activity in primary motor cortex is elevated in PD (Crowell et al., 2012), as well as excessively coupled to the phase of low-frequency cortical rhythms (phase-amplitude coupling) (de Hemptinne et al., 2013). Thus, analysis of the interaction of basal ganglia spiking with cortical broadband could yield great insight into basal ganglia-cortex dynamics in PD.
To address these questions, we recorded primary motor cortex (M1) arm area LFPs simultaneously with STN unit discharge, in awake patients undergoing deep brain stimulator (DBS) implantation surgery. M1 LFPs were recorded using subdural electrocorticography (ECoG), a technique that combines excellent spatial and temporal resolution with sufficient signal amplitude to resolve high-frequency broadband activity, in addition to low-frequency rhythms, relatively free of artifact. To provide a comparison group, we performed similar recordings in patients with primary craniocervical dystonia who had minimal symptoms in the contralateral arm. We tested the hypothesis that in PD, STN single-unit discharge at rest is excessively synchronized to cortical oscillations, to cortical broadband gamma activity, and to epochs of cortical phase-amplitude coupling. Our findings support a model of the basal ganglia-thalamocortical loop in PD in which cortical 4–30 Hz rhythms are phase-locked to gamma activity, and in which phase modulated “waves” of gamma activity in motor cortex drive STN hyperactivity.
Materials and Methods
Subject recruitment and clinical characterization.
Study subjects were recruited from a population of movement disorders patients scheduled to undergo DBS implantation at one of two campuses: the University of California at San Francisco (UCSF), or the San Francisco Veterans Affairs Medical Center (SFVAMC). Subjects had a diagnosis of idiopathic PD, or primary cervical or craniocervical dystonia, confirmed by a movement disorders neurologist (J.L.O. or N.B.G.). The dystonia patients included in this study were simultaneously participating in a clinical trial to evaluate the effectiveness of STN DBS in dystonia (Ostrem et al., 2011). Informed consent was obtained before surgery under a protocol approved by the UCSF/SFVAMC Institutional Review Board, according to the Declaration of Helsinki. During the consent process, it was explained to each prospective subject that the temporary intraoperative placement of the cortical recording strip was performed solely for research purposes. Potential study subjects were characterized within 60 d before surgery using the following standardized rating scales: for PD patients, the Unified Parkinson's Disease Rating Scale part III (UPDRS-III) after withdrawal of anti-parkinsonian medications for 12 h; for dystonia patients, the Burke–Fahn–Marsden dystonia rating scale (BFMDRS) movement score, and the Toronto Western Spasmodic Torticollis Rating Scale (TWSTRS) severity score. For PD patients, preoperative UPDRS-III subscores for the hemibody contralateral to brain recordings were used in some analyses for correlations with physiological data. However, tremor was also scored for each individual recording as “present” or “absent” intraoperatively based on off-line visual inspection of limb electromyography (EMG) and accelerometry (detailed further below). Preoperative medications in PD were summarized as levodopa-equivalent doses (LED) using the following conversion factors: ropinirole × 20; pramipexole × 100; levodopa with decarboxylase inhibitor × 1; controlled release levodopa with decarboxylase inhibitor × 0.7; levodopa with decarboxylase and COMT inhibitor × 1.3 (Wenzelburger et al., 2002).
Subject inclusion criteria were as follows: age 21–75 years, normal brain magnetic resonance imaging (MRI) examination, sufficient disease severity in the setting of optimal medical management to justify treatment by DBS, and ability to cooperate during awake neurosurgery. Further criteria for the primary dystonia group: dystonia affecting predominantly cervical or craniocervical muscles with minimal arm involvement, and no treatment with injected Botulinum toxin for 3 months before surgery. Twenty-nine PD patients and six dystonia patients met study criteria and were included.
Placement of subdural ECoG electrodes and STN microelectrodes.
All subjects underwent typical procedures for planning and stereotactic surgical placement of deep brain stimulator electrodes into the STN (Starr et al., 2002). On the surgical planning MRI, the STN target was identified as a signal hypointensity on T2-weighted fast spin-echo imaging lateral to the anterior border of the red nucleus, typically 11–13 mm from the midline. For targeting M1 for placement of the ECoG array, we identified a point on M1, 3 cm from the midline, based on anatomic identification of the central sulcus (Fig. 1A). The intention was to target the arm-related area of M1, slightly medial to the “hand knob” (Yousry et al., 1997), in the same parasagittal plane as the typical surgical entry site for placement of STN DBS electrodes. A radio-opaque marker was stereotactically placed on the scalp over this location. When bilateral DBS implantations were planned, cortical recording was performed on one side only, typically the side with the clearest anatomic demarcation of the central sulcus. After drilling of frontal burr holes, and opening of the dura mater, we placed a six-contact subdural ECoG strip (3 mm contacts, 1 cm spacing) (Ad-Tech or Integra) on the brain surface and directed it posteriorly in a parasagittal plane under fluoroscopic control to provide contacts covering both M1 and primary sensory cortex (S1). The locations of the burr holes were determined solely by the selection of the safest entry point for the intended DBS trajectory, and no additional skull or scalp exposure was needed for ECoG strip placement. After placement of the ECoG strip and guide tube for microelectrode recording, the burr hole was sealed with a fibrin sealant.
Method of localization of ECoG electrode contacts and of STN recording site. A, Intended location (white arrow) of the center of the ECoG recording strip, immediately anterior to the central sulcus (CS), on preoperative T1-weighted planning MRI, axial view. B, Actual location of the recording contacts on the ECoG strip, from intraoperative CT computationally fused with the preoperative planning MRI, on parasagittal view 3 cm from the midline. C5 is the contact best localized to M1 by MRI in this example. C, Somatosensory evoked potential (from stimulation of the median nerve) for three adjacent contact pairs. Reversal of the N20 potential at contact pairs more posterior than C5–C6 confirms C5 as the contact best localized to M1. D, Location of the STN DBS electrodes. The intraoperative CT, showing the lead artifacts (white arrows), is computationally fused with the preoperative axial T2-weighted MRI, which shows the STN as a region of T2 hypointensity.
Confirmation of ECoG electrode location.
ECoG contact locations were confirmed anatomically using intraoperative CT (iCT) (O-arm, Medtronic) (Shahlaie et al., 2011) or lateral fluoroscopy. iCT images were computationally fused to the preoperative planning MRI using standard surgical planning software (Framelink 5.1, Medtronic), and windowed so as to visualize each ECoG contact with respect to underlying MRI anatomy (Fig. 1B). For cases where lateral fluoroscopy only was used, a lateral x-ray documented the location of the ECoG contact in the anterior–posterior direction with respect to the radio-opaque scalp marker.
For physiological confirmation of electrode location, we recorded somatosensory evoked potentials (SSEPs) generated by median nerve stimulation. The stimulation parameters were as follows: pulse frequency, 2 Hz; pulse width, 200 μs; pulse train length, 160; current 25–35 mAmp (typically 50% higher than the threshold for visible thumb twitch). Signals recorded from the ECoG strip during the SSEP were sampled at 5000 Hz, bandpass filtered (1–500 Hz), and amplified × 7000. Signal averaged SSEPs for each adjacent contact pair (computationally re-montaged) were then visually inspected to determine the M1 location by reversal of the N20 waveform. For each bipolar pair, the more posterior contact was used as the “active” electrode, whereas the more anterior was the “reference”. The posterior contact of the most posterior bipolar contact pair that showed a negative going waveform was considered to be the contact best localized to M1 (Fig. 1C). There was high concordance between anatomic and physiologic determinations of contact position with respect to the central sulcus.
Cortical LFP and STN single-unit recording.
Cortical LFPs and STN units were recorded in a state of alert rest, using the Guideline 4000 system (FHC) or the Alpha Omega Microguide Pro (Alpha Omega). These are customized 14-channel data collection systems approved by the United States Food and Drug Administration for human use. Five-channel cortical LFP recordings were obtained from each of the five most posterior ECoG contacts (C1–C5), each differentially referenced to the most anterior contact (C6). A needle electrode in the scalp served as the ground. LFP signals were bandpass filtered 1–500 Hz, amplified × 7000, and sampled at 1000 Hz (Guideline 4000 system) or 1500–3000 Hz (Alpha Omega). A 60 Hz notch filter was used to reduce line artifact. Data collected with the Alpha Omega system were downsampled to 1000 Hz (MATLAB function resample). All subjects had been sedated with propofol for the initial surgical exposure, but data collection was performed at least 30 min after stopping propofol. This is sufficient time for all neuronal effects of this agent to be eliminated (Fechner et al., 2004; Raz et al., 2008). Simultaneous EMG (biceps brachii, extensor carpi radialis, and gastrocnemius) and arm accelerometry were performed to confirm the absence of voluntary movement, and presence or absence of tremor, during recording. All raw, unfiltered EMG and accelerometry recordings were visually inspected off-line, and if 4–6 Hz oscillations were present in at least one EMG or accelerometry trace for a minimum continuous duration of 10 s, the corresponding single-unit/M1 LFP recording was scored as “tremor present”; otherwise it was scored as “tremor absent”. When present, 4–6 Hz oscillations typically were seen during most or all of the recording, and scoring was based on concurrence between two examiners. Because tremor can occur transiently, some recordings from a single subject could be scored as tremor present, whereas others were tremor absent.
STN neuronal recordings were obtained using glass-coated platinum/iridium microelectrodes with impedance 0.4–1.0 MΩ at 1000 Hz (FHC or Alpha Omega) as previously described (Starr et al., 2002). Microelectrodes were advanced into the brain using a motorized or manual microdrive (Alpha Omega or Elekta). In a typical surgical case, 1–2 microelectrode penetrations were made serially through the STN on each side, separated by 2–3 mm. Signals were bandpass filtered (300 Hz to 3 kHz), amplified, played on an audio monitor, displayed on an oscilloscope, and digitized (20 kHz sampling rate). Cells were recorded at approximately every 300–800 μm along each trajectory, and each recording lasted 30–120 s. Neurons were screened for movement-related activity based on audible changes in action potential discharge evoked by passive (investigator-initiated) contralateral limb movements (i.e., shoulder, elbow, wrist, hip, knee, and ankle joints). Proprioceptive responsiveness of a neuron in relation to movements of one or more joints was determined by concurrence between the examiner, and one or more operating room staff based on audiovisual assessments of the response. Only neurons recorded in the dorsal 4 mm of the nucleus, in a region where movement related activity was detected, and in which a minimum of 1000 spikes were digitized, were used for subsequent analysis.
Confirmation of STN localization.
Most STN units (89%) were recorded on the same trajectory as the final DBS electrode placement. Postoperative MRI, or CT computationally fused to the preoperative MRI (Framelink 5.1 software, Medtronic) (Fig. 1D), was performed to verify appropriate localization of the DBS electrode in the STN.
Single-unit discharge properties.
Digitized spike trains were imported into off-line spike-sorting software (Plexon) for discrimination of single populations of action potentials by principal components analysis (Adamos et al., 2008). The software generated a record of spike times (subsequently reduced to millisecond accuracy) for each action potential waveform detected. The first 1000 spike times were used to calculate discharge rate and oscillatory activity in the 0–200 Hz range. Analyses were performed using MATLAB software (MathWorks). Neuronal data were accepted as single-unit data in this study only if action potentials could be discriminated with a high degree of certainty, as indicated by the presence of a clear refractory period of at least 3 ms in the interspike interval (ISI) histogram. Neurons whose action potential morphology varied greatly in synchrony with the cardiac cycle were excluded. Oscillations in the spike train at 0–50 Hz were evaluated using the “global spike shuffling” method (Rivlin-Etzion et al., 2006) to eliminate the artifactual autocorrelations that arise from the neuronal refractory period. Spike time stamps were converted to a data stream consisting of 1 ms bins in which the occurrence or absence of a spike was represented by 1 or 0, respectively, in that bin. A 2048-point fast Fourier transform with Hanning window was used, resulting in a spectral resolution of 0.5 Hz. Similar analysis was performed on “control” data in which ISIs had been randomly shuffled 100 times (Rivlin-Etzion et al., 2006). Statistically significant peaks in the spike train data (after normalization with the spectrum of the shuffled data) were determined by using the 300–500 Hz part of the spectrum as the control segment and its SD was used as a measure of random fluctuations in the spectrum. Each frequency point between 0.25 and 200 Hz was then checked for deviation from the expected power, at a significance level of p = 0.0025, after correction for multiple (100) comparisons. Only 4–50 Hz oscillations were plotted as few occurred outside of this range.
LFP power spectra.
For analysis of cortical LFPs, the raw strip recordings were computationally re-referenced to adjacent bipolar pairs by subtraction of the common reference. The M1 contact was referenced to either the contact immediately anterior or the contact immediately posterior to it, depending on which pair showed a larger spike-timed average waveform (described further below). The M1 LFP power spectrum was calculated using the Welch method (MATLAB function pwelch) with a 512-point FFT, for 1.96 Hz resolution, using the entire LFP recording (duration of recordings 34.3–140.0 s).
Data analysis: STN spike synchronization to M1 LFP oscillations.
The first 1000 consecutive spike times of each recording were used for all analyses of spike-cortex synchronization. To provide a sensitive measure of spike-LFP interactions and to calculate latencies between the timing of spikes and peaks in cortical oscillatory activity, we analyzed the spike-timed average (STA) of the M1 LFP. The STA was computed in a 1 s window centered on the time of occurrence of STN spikes. We used the STA to calculate an STN spike-M1 LFP oscillation modulation index MIO as follows (see Fig. 3B–D). (1) The peak frequency of the STA was determined from the maximum of its power spectrum (512-point FFT, MATLAB function pwelch, 1.96 Hz resolution) (see Fig. 3B, inset). (2) We bandpass filtered the STA with a frequency band of width 4 Hz, centered on the peak frequency [100 tap FIR filter implemented using MATLAB function filtfilt for bidirectional filtering to avoid filter-induced phase shifts (Moran et al., 2008)]. (3) We calculated the amplitude envelope of the bandpass-filtered STA time series as the magnitude of its (complex valued) analytic signal (see Fig. 3C). The analytic signal was generated from the sum of the bandpass-filtered STA and its Hilbert transform (MATLAB function hilbert). (4) We restricted the analysis to a 400 ms window centered on time 0. The maximum value of the STA amplitude envelope was determined, and a mean value (Mtest) was computed using a 100 ms window around the peak—this was considered to be a “raw” measure of spike-M1 rhythm synchronization. For the same time period, a surrogate mean amplitude envelope (Mrev) and its SD (SDrev) were computed from the same cortical LFP data, but using temporally reversed spike times (see Fig. 3D). Mrev and SDrev were computed from a single time-reversed dataset (400 samples over a 400 ms window) providing a simple measure of STA fluctuations in the absence of synchronization. We defined the STN spike-M1 oscillation modulation index MIO as (Mtest − Mrev)/SDrev, thus generating a z-scored index of the magnitude of synchronization between simultaneously recorded STN spikes and M1 LFPs. (5) Some STAs contained two frequency components, as indicated by two peaks in their power spectra; in these cases, an MIO was also defined for the secondary frequency.
MIO values >3 were considered to be statistically significant. This threshold resulted in the occurrence of false-positives (recordings in which the reverse time STA analytic amplitude Mrev also crossed the 3.0 SD level at some time within the 400 ms STA test window) at a frequency of <5%. For all STAs that exceeded the threshold of significance, we tabulated: (1) the frequency of peak spectral power, calculated as described in step 1 above (see Fig. 3B); (2) the instantaneous phase value of the bandpass-filtered STA at spike time (time 0) (see Fig. 3C); (3) the latency from time 0 to the time of the peak of the amplitude envelope of the bandpass-filtered STA (see Fig. 3C); and (4) the duration of the interval over which the amplitude envelope of the STA crossed and stayed above the threshold of 3 SDs (see Fig. 3D). Of note, the exact time at which an STA amplitude envelope crosses its significance threshold is in part related to the frequency of the underlying oscillation, not just the strength of the STN-M1 interaction. Therefore, analyses of latency and duration of spike-LFP synchronization were only used as general descriptors, not for detailed comparison between subject groups.
To confirm the overall frequency characteristics of spike-M1 synchronization obtained by analyses of STN spike-timed averages of the M1 LFP, we also calculated STN spike-cortical LFP coherence using traditional methods for determining coherence between a point process and a field potential, described by Halliday et al. (1995).
Data analysis: STN spike synchronization to M1 broadband gamma activity.
Because low frequencies (<50 Hz) predominate in the LFP power spectrum and its STA, we used additional methods to study synchronization between STN spikes and cortical broadband gamma activity. First, we visualized spike-broadband gamma interactions by plotting the spike-timed scalogram of the M1 LFP, in 1 s windows centered on the time of occurrence of STN spikes, from 0 to 300 Hz, using a wavelet decomposition (see Fig. 4A). A Morlet wavelet of 5 cycles in width was convolved with the M1 LFP time series to estimate an amplitude and phase of the signal at each frequency for every point in time.
Next, we calculated an STN spike-M1 gamma modulation index MIgamma using methods similar to those above for MIO, but first isolating the 50–200 Hz portion of the cortical LFP spectrum (see also Fig. 4B–E). Before signal averaging, we bandpass filtered the raw cortical LFP at 50–200 Hz to isolate broadband gamma (100 tap FIR filter, implemented using MATLAB function filtfilt), and calculated the amplitude envelope of broadband gamma time series, from the magnitude of its analytic amplitude. This time series is referred to as gamma_AE. We computed the spike-timed average (STA) of gamma_AE, in 1 s windows centered on the time of occurrence of STN spikes (N = 1000 spikes). This STA waveform is indicative of the degree to which cortical broadband gamma activity is synchronized to the timing of STN spikes. We noted that, when spike-synchronized peaks in broadband gamma occurred, the amplitude of cortical broadband gamma was invariably modulated by the phase of a low-frequency rhythm between 4 and 30 Hz (see Fig. 4B). We therefore defined the spike-gamma modulation index MIgamma in a manner that accounted for this behavior, using the amplitude envelope of the STA of gamma_AE. The method for determining MIgamma was identical to steps 1–5 described above for calculating the spike-M1 oscillation modulation index MIO, except that the analysis used the STA of gamma_AE rather than the STA of the unfiltered M1 LFP.
MIgamma values > 3 were considered to be statistically significant. For all spike-LFP pairs with a significant MIgamma, we calculated: (1) the frequency of peak spectral power of the STA of gamma_AE (see Fig. 4B); (2) the instantaneous phase value of the bandpass-filtered STA of gamma_AE at spike time (time 0) (see Fig. 4D); (3) the latency from time 0 to the time of the peak of the amplitude envelope of the bandpass-filtered STA of gamma_AE (see Fig. 4D); and (4) the duration of the interval over which the amplitude envelope of the STA of gamma_AE crossed and stayed above the threshold of 3 SDs (see Fig. 4E).
For the subset of recordings that had modulation indices >3 and that had >2000 consecutive spikes available for analysis, we checked the effect of spike number (500, 1000, or 2000) on the modulation indices MIO and MIgamma. For most recordings, spike-M1 synchronization reached statistical significance using a spike number of 1000, justifying the use of trains of 1000 spikes for most analyses in the study.
Statistical analysis of grouped data.
Parameters calculated from the STN spike-cortical LFP synchronization analyses were summarized by their means and SDs, or medians and ranges depending on the normality of the data (tested using MATLAB function lillietest). Of note, for most study subjects, multiple STN-M1 paired recordings were sampled, from different locations within the motor territory of the STN (Tables 1 and 2). We did not assume a priori that different recordings from the same subject were statistically independent. Therefore, for comparisons between disease groups (see Fig. 5A), results from standard statistical tests (Wilcoxon rank sum test) were also checked using a general estimation equations method (GEE, SPSS Software), with subject number as a within-subjects variable to account for potential clustering of within-subjects data (Hanley et al., 2003). Similarly, correlations of physiological measures with symptom severity (see Figs. 3E, 4F) were performed both with linear regression (assuming independence of all samples), and using a GEE model that did not assume sample independence within subjects. Phase distributions were analyzed for their difference from a uniform distribution using the Rayleigh test (Moran et al., 2008). Categorical data were analyzed using χ2 or Fisher's exact tests.
Characteristics of Parkinson's disease subjects
Characteristics of primary dystonia patients
Results
Study subjects
Characteristics of the study subjects, including preoperative medications, are provided in Tables 1 and 2. There were 29 subjects in the PD group, and six in the dystonia group. Primary dystonia patients had mainly craniocervical involvement, with minimal symptoms in the arm contralateral to the side of recording in most subjects. Two dystonia subjects had action-induced arm tremor or writer's cramp, resulting in nonzero BFMDRS scores for the contralateral arm (Table 2), but did not have tremor or dystonic posturing of limbs at rest.
Oscillations in STN unit discharge and M1 LFPs in PD
In PD patients, a total of 116 simultaneous recordings of STN single units with M1 LFPs were analyzed. Forty-six of these recordings occurred during contralateral arm tremor, whereas 70 occurred during no detectable contralateral arm tremor. The spike firing rate (mean ± STD) was 35.6 ± 14.1 Hz. Twenty-eight units had significant oscillations in the spike train in the 4–50 Hz range (Fig. 2A,C). For subsequent analyses, oscillations were further subcategorized as beta frequency (12–30 Hz), theta-alpha frequency (4–12 Hz, which includes tremor frequency and its first harmonic), or neither (Table 3). The distribution of single-unit oscillation frequencies was bimodal, with one beta cluster and one theta-alpha cluster, consistent with prior reports (Levy et al., 2000; Rodriguez-Oroz et al., 2001; Moran et al., 2008). Two units oscillated in both frequency bands (Fig. 2A, middle example). Oscillations near tremor frequency were found in some recordings in the absence of detectable tremor during the recording, but were more likely to occur when tremor was present (details in Table 3).
Oscillatory activity in STN single-unit discharge and in the M1 LFP between 4 and 50 Hz, in PD. A, Examples of units with tremor frequency, beta frequency, and both types of oscillatory activity. A 1 s recording of neuronal activity is shown with the frequency spectrum for the spike train below. The gray dotted line is the significance threshold, and gray circles indicate statistically significant peaks in the power spectrum (Rivlin-Etzion et al., 2006). B, Simultaneously recorded M1 LFPs, and their power spectra are shown in bottom row for each example. Vertical arrows show peaks in the power spectrum in theta or beta ranges. C, D, Distribution of single-unit oscillation frequencies (C) and peak frequencies in the M1 LFP power spectrum (D) across all PD subjects, segregated according to the presence or absence of tremor during the recording.
Proportion of recordings with different types of oscillatory behavior, categorized by presence of or absence of tremor
Power spectra of the M1 LFP (Fig. 2B) were inspected for peaks (local maxima) in the beta and tremor bands. All but one M1 recording had a spectral peak in the beta band, whereas 68 also had a separate peak in the theta or alpha frequency bands (Fig. 2B,D). The presence or absence of contralateral limb tremor during the recording was not predictive of the presence or absence of a tremor frequency peak in the M1 LFP power spectrum (Fig. 2D; Table 3).
Synchronization of STN spikes to motor cortex oscillations at 4–50 Hz in PD
Synchronization of STN single units to cortical rhythms was examined primarily by analysis of STN spike-timed averages of the cortical LFP. The STA waveform is indicative of the degree to which cortical oscillations are phase synchronized to the timing of STN spikes. An illustrative example is shown in Figure 3A,B. When compared across all five cortical LFP recording pairs (Fig. 3B), spike-cortex synchronization was always maximal in a contact pair that covered M1. The methods of quantifying the frequency, phase, and magnitude (modulation index MIO), of the STN spike synchronization to M1 LFP oscillations, are shown in Figure 3B–D. The magnitude of spike-M1 LFP synchronization correlated with the severity of contralateral bradykinesia (using preoperative UPDRS-III bradykinesia scores) (Fig. 3E), and this was true whether or not within-subject recordings were treated as statistically independent (p value for linear correlation of independent samples of 0.015 with Pearson's r of 0.224, versus p value from the GEE model of 0.001 with a beta coefficient of 0.37).
Synchronization of STN spiking with cortical LFP oscillations in PD. A, Example of simultaneously recorded STN unit activity (top trace) and M1 LFP (bottom trace), >1 s. B, Example of STA of the cortical LFP, for all five contact pairs from most posterior (top trace) to most anterior (bottom trace; N = 1000 spikes). The inset shows the power spectrum of the M1 STA. C, Method of determining phase and latency of M1 STA with respect to STN spike time. Black line shows bandpass-filtered STA, red line shows its amplitude envelope, and vertical dashed line marks spike time. Open circle, the zero crossing of the bandpass-filtered STA (phase = 1.4 radians). Closed circle, peak of the amplitude envelope, occurring at a latency of 7 ms before spike occurrence. D, Method of determining the spike-M1 LFP oscillation modulation index MIO. Red line is the amplitude envelope of the bandpass-filtered STA, with dotted horizontal bar showing average value (during a 100 ms window centered on its peak). Blue line is the amplitude envelope of the reverse time STA, with horizontal lines showing mean ± 3 SDs during a 400 ms window. The M1O for this unit-M1 pair is 9.2 in z-score units, and the STA frequency is 23.4 Hz. E, Positive correlation between MIO and preoperative UPDRS bradykinesia scores for the contralateral body. F, Box plots showing the distribution of time of onset, and time of peak effect, of spike-M1 gamma synchronization. G, Histogram showing the distribution of peak frequencies for the STN spike time averaged M1 LFP in PD patients with and without tremor during the recording. H, Phase histogram, showing the distribution of phases for the spike time averaged M1 LFP waveform with respect to STN spike time. Data in G and H are only for recordings with MIOs that reached statistical significance. I, Bar graph showing the numbers of recordings with and without significant MIO, segregated by the responsiveness of STN cells to contralateral limb movement. NR, No response.
For those spike-M1 pairs with significant modulation (defined as MIO > 3, 56 recordings), the latency, frequency, and phase relationships between spike and LFP are shown in Figure 3F–H. Three recordings had significant spike-LFP modulation at two different frequencies. The median onset of spike-LFP synchrony was 103 ms before the spike time (p < 0.001 for difference from zero, nonparametric sign test) and the median time of peak modulation (defined as shown in Fig. 3C) was 18.5 ms before the spike time (not significantly different from zero by the sign test) (Fig. 3F). The distribution of STN-M1 LFP synchronization frequencies, derived from the power spectra of the spike-timed averages, showed clusters in the theta and beta bands. (Fig. 3G). STN-M1 LFP synchronization at or near tremor frequency (4–12 Hz) was much more prominent when limb tremor was present during the recording (Table 3). Similar results were found when STN spike–M1 LFP coherence was analyzed by the Halliday method rather than by using spike time averages of M1 (data not shown). The distribution of phases for STN-M1 synchronization was significantly different from a uniform distribution (p = 0.001, Rayleigh test), with a mean phase angle of 1.6 radians, corresponding to the falling slope of the bandpass-filtered M1 LFP sinusoid (Fig. 3C,H). STN spike synchronization to the M1 LFP was not restricted to neurons that had prominent arm related activity during intraoperative somatosensory examination, but was also found for neurons whose discharge was primarily modulated by leg movement or had no detectable response to passive movement (p = 0.83, χ2) (Fig. 3I), notwithstanding the fact that cortical recordings were made exclusively over the arm territory of M1.
STN spike-cortical gamma synchronization
The use of the spike-timed LFP average for analysis of spike-M1 synchrony (Fig. 3) mainly explores the low-frequency range of the LFP power spectrum (<50 Hz), as these frequencies predominate in the LFP and in its STA. However, we were also interested in synchronization between STN spikes and cortical broadband gamma activity (50–200 Hz), as the latter is thought to reflect underlying cortical spiking activity (Manning et al., 2009), and therefore offers a means of exploring the relationship between STN spiking and cortical spiking. To visualize STN spike-M1 gamma interactions, we plotted spike-timed scalograms, illustrated for one representative example in Figure 4A. In the scalogram, spectral power is normalized separately at each frequency (to the mean of its wavelet convolution), allowing visual examination of gamma frequencies that are too small to observe in the spike-timed LFP average illustrated in Figure 3.
Synchronization of STN spiking with M1 broadband gamma activity in PD. A, Example of spike-timed scalogram, showing cortical activity aligned to the time of STN spikes. Color-scale represents the amplitude of oscillatory activity, normalized separately for each frequency. B, One-second segment of M1 ECoG, bandpass filtered to isolate broadband gamma component (50–200 Hz) (blue line), with its amplitude envelope (red line). C, Spike-timed average of the M1 broadband amplitude envelope (1000 spikes), showing amplitude modulation of the gamma component (phase-amplitude coupling) preceding and overlapping with the timing of STN spikes. The inset shows the PSD of the spike-timed average of the M1 gamma amplitude. D, Method of determining the latency and phase of the M1 gamma amplitude envelope with respect to STN spike time. Black line shows the STN spike time average of the broadband gamma component (50–200 Hz) of the M1 LFP (after bandpass filtering around its peak frequency). Open circle marks the phase of the signal at time 0 (1.7 radians in this example). Red line shows the amplitude envelope of the spike time average of the broadband gamma component. Filled circle shows the peak of amplitude envelope, occurring at a latency of 58 ms before spike occurrence in this example. E, Method of determining the spike-gamma modulation index MIgamma. Solid red line is the solid red line in D, with horizontal dotted red bar showing average value over a 100 ms interval around the peak. Blue line is the amplitude envelope of the reverse-timed spike time average of M1 broadband gamma activity, with horizontal dotted blue lines showing mean ± 3 SDs (MIgamma is 8.0 in this example). F, Negative correlation between MIgamma and preoperative UPDRS resting tremor scores for the contralateral body. G, Box plots showing the distribution of time of onset, and time of peak effect, of spike-M1 gamma synchronization. H, I, Histograms showing the distribution of frequencies (H) and phases (I) for the oscillatory amplitude envelope of the spike time average of the M1 gamma component, for recordings whose MIOs reached statistical significance. Data in H are segregated by the presence or absence of tremor during the recording. J, Bar graph showing the number of recordings with and without significant MIgamma, segregated by the responsiveness of STN cells to contralateral limb movement. NR, No response.
Inspection of spike-timed M1 scalograms revealed synchronization of broadband gamma with STN spikes, in a highly distinctive pattern: waves of increased M1 broadband gamma activity, alternating with troughs of broadband gamma. The occurrence of this pattern in the scalogram implies that there are epochs of cortical phase-amplitude coupling whose phase frequency is itself synchronized to the timing of STN spikes. The highest amplitude gamma waves typically preceded the time of occurrence of STN spikes. This oscillating pattern of broadband gamma activity had an amplitude envelope typically in the beta, alpha, or theta range. The frequency range of broadband gamma that was most prominently synchronized to STN spiking was 50–200 Hz (Fig. 4A), so this region of the spectrum was used for further quantitative analysis.
Figure 4B–E shows the method of quantifying the parameters of STN spike–M1 gamma synchronization, including the frequency and phase of the amplitude envelope for gamma activity, the latency from peak gamma modulation to STN spike time, and the overall strength of the spike-gamma interaction (modulation index MIgamma). The magnitude of spike-M1 gamma synchronization did not correlate with lateralized total preoperative UPDRS scores or with contralateral bradykinesia scores, but did have a modest negative correlation with contralateral tremor scores (Fig. 4F). These results were invariant to the statistical method used for correlation (treating all recordings as independent samples, p = 0.038 with r of −0.193, or treating within-subjects recordings as nonindependent samples in a GEE model, p = 0.043 and beta coefficient of −0.226, for the inverse correlation with tremor). For those spike-M1 pairs with significant modulation of M1 gamma (defined as MIgamma > 3, 46 recordings), the latency, frequency, and phase relationships for synchronization of STN spikes to the M1 gamma amplitude envelope are shown in Figure 4G–I. The median onset of spike-gamma synchrony was 107.1 ms before the spike time (p < 0.001 for difference from zero, nonparametric sign test) and the median time of peak modulation (defined as shown in Fig. 4D) was 48 ms before the spike time (p = 0.05 for difference from zero, sign test, Fig. 4G). The distribution of frequencies for the spike-synchronized M1 gamma envelope, derived from the power spectra of the spike-timed average of gamma activity, shows that the gamma amplitude envelope oscillated most prominently in the theta and beta bands (Fig. 4H). Spike-gamma envelope synchronization at any frequency, and especially at beta frequencies, was more likely to occur in the absence of contralateral tremor (Table 3). This contrasts with the relationship between spike-LFP synchronization and contralateral tremor (Table 3), and suggests a relationship between symptom profile and the type of synchronization (MIO vs MIgamma), not just the frequency, of basal-ganglia-cortex interactions. The distribution of phases for STN-M1 gamma envelope synchronization showed variable phase relationships that did not differ from a uniform distribution (p > 0.05, Rayleigh test) (Fig. 4I). STN spike synchronization to M1 gamma was not restricted to neurons that had prominent arm related activity during intraoperative somatosensory examination (p = 0.72, χ2) (Fig. 3J).
Disease-specific patterns of STN-M1 synchronization
Although normal control data for invasive recordings cannot be obtained from humans, we were able to assess the disease specificity of our findings in PD subjects by comparison with patients with primary craniocervical dystonia (Fig. 5), who also underwent placement of stimulators into the STN (Ostrem et al., 2011). Twenty-six spike-LFP pairs from six dystonia patients were available for analysis. For STN synchronization to M1 oscillations in the 4–50 Hz range, the magnitude of synchronization (MIO) was less than in PD (p = 0.005 for Wilcoxon rank sum test and 0.009 for GEE model) (Fig. 5A, left). Six spike-LFP pairs in dystonia patients had significant synchronization (MIO > 3), in the theta, alpha, or beta ranges (Fig. 5B, left), but the proportion of all recordings with significant beta range spike-LFP synchronization was less than for PD (4 of 26 in dystonia vs 41 of 116 in PD, p = 0.05, χ2). For STN synchronization to M1 gamma activity, six paired recordings in dystonia also showed waves of amplitude modulated broadband activity, slightly preceding the time of occurrence of STN spikes, one of which occurred at two different frequencies of the sinusoidal amplitude envelope (Fig. 5B, right). However, the pattern of spike-gamma synchronization was strikingly different from that seen in PD, in that the frequency of the sinusoidal amplitude envelope for broadband gamma was never in the beta band in dystonia patients but was commonly in the beta band in PD patients (0 of 26 in dystonia vs 23 of 116 in PD, p = 0.008 by Fisher's exact test).
Comparison of STN spike-M1 synchronization in PD versus primary craniocervical dystonia. A, Box plots summarizing MIO (left) and MIgamma (right) over all STN-M1 recordings in PD versus dystonia. p values are from the Wilcoxon rank sum test. B, Histograms showing the distribution of frequencies for STN spike-M1 LFP synchronization (left) and for STN spike-M1 gamma synchronization (gamma amplitude envelope, right). Statistical analysis of the proportion of recordings with significant synchronization is provided in the results section. Data in B are only for recordings with an MIO or MIgamma that reached statistical significance (>3 in z-score units).
Relation between spike-M1 LFP rhythm modulation and spike-M1 gamma modulation
Two forms of spike-M1 synchronization were studied here: spike synchronization to M1 rhythms (4–50 Hz), and spike synchronization to M1 gamma (50–200 Hz). Several lines of evidence suggest that these represent distinct phenomena. First, across all recordings with a significant MIO or significant MIgamma, the value of these modulation indices were not correlated (Fig. 6A,B). Second, although the spike-gamma modulation had an oscillating amplitude envelope (Fig. 4A–C), the frequency of the amplitude envelope for the STN spike-M1 gamma modulation did not correlate with the frequency for spike-M1 rhythm modulation (Fig. 6B). Thus, the mechanisms for generating these two modes of synchronization may be at least partly independent.
Relationship between STN spike-M1 LFP synchronization and STN-M1 gamma amplitude envelope synchronization in PD. A, Lack of correlation between the magnitudes of MIO and MIgamma, for all recordings with MIO or MIgamma (or both) >3. B, The amplitudes (top plot) and frequencies (bottom plot) for the two types of synchronization (MIO vs MIgamma) are plotted for a subset of five PD subjects who each had at least 4 different units showing significant synchronization of at least one type. Data from each subject are color coded, to show that even within a given subject the oscillation frequencies for individual spike-LFP pairs were variable. In this subset, as in A, there was no correlation between MIO and MIgamma with respect to either amplitude (p = 0.31) or preferred frequency (p = 0.32). In the bottom plot, the identity line is drawn to facilitate comparison between preferred frequencies for the different types of synchronization.
Discussion
To examine cortex-basal ganglia interaction in humans with movement disorders, we recorded M1 local field potentials simultaneously with STN single-unit activity in patients undergoing STN DBS implantation for Parkinson's disease, as well as in a comparison group with primary craniocervical dystonia. LFPs represent summed, synchronized subthreshold activity in presynaptic and postsynaptic elements near the recording electrode, and are increasingly used in studies of neuronal synchronization in the motor system in normal, as well as disease states (Murthy and Fetz, 1996; Goldberg et al., 2004; Hammond et al., 2007). In general, the phase of LFP oscillations biases the probability of spike firing. In the normal state, this modulation is dynamic and task specific, representing a flexible mechanism whereby LFP oscillations can link functionally related brain areas (Fries, 2005; Leventhal et al., 2012).
STN spike-M1 rhythm synchronization
Excessive neuronal oscillatory activity within specific basal ganglia nuclei is a well documented feature in parkinsonian animal models as well as in human PD (Dostrovsky and Bergman, 2004; Hammond et al., 2007). Similarly, simultaneously recorded neurons within the subthalamic nucleus or internal globus pallidus have oscillations that are strongly synchronized (Nini et al., 1995; Levy et al., 2002a; Hanson et al., 2012), and basal ganglia neurons tend to be synchronized to LFPs recorded in the same structure (Kuhn et al., 2005; Moran et al., 2008). However, the relation of basal ganglia oscillations to cortical activity is less well studied. In humans, interactions between cortex and basal ganglia at fast time scales have been primarily studied by recording basal ganglia LFPs in conjunction with scalp electroencephalography (Williams et al., 2002; Fogelson et al., 2006; Lalo et al., 2008; Eusebio et al., 2009) or magnetoencephalography (Litvak et al., 2011). These studies have elucidated specific patterns of cortex-basal ganglia coherence, but in general have not used nonparkinsonian comparison groups.
Here, we provide the first demonstration of the synchronization of basal ganglia spiking with cortical oscillations in humans with PD. This behavior is relatively specific to the parkinsonian state, as it is less prominent in a different movement disorder of basal ganglia origin, primary craniocervical dystonia. We show that synchronization occurs across a broad range of low-frequency rhythms (theta, alpha, and beta). Spike-LFP synchronization at tremor (theta) frequency is much more likely when tremor is present during the recordings (Fig. 3G). However, tremor frequency oscillations in spike discharge and M1 LFPs do occur in the absence of tremor (Fig. 2), suggesting that neuronal oscillatory phenomena at tremor frequency are not exclusively a consequence of tremor-frequency peripheral input. It is possible that tremor frequency discharges at different points in the circuit represent a “parkinsonian endophenotype,” but that the clinical manifestation of tremor specifically requires synchronization of these oscillations between basal ganglia and cortex.
Our findings support the evolving view that excessive synchronization between basal ganglia and cortex is a prominent feature of PD. We provide evidence that the parkinsonian state alters the statistical relationship between LFP phase and spike timing in global brain networks, not just in local basal ganglia or cortical circuits (Goldberg et al., 2004; Gatev and Wichmann, 2009; Li et al., 2012). Because many brain functions appear to depend on this mechanism for task performance (Fries, 2005), its disruption in a critical motor circuit could result in widespread impairment of motor functions. Our findings in humans are broadly consistent with several rodent studies of parkinsonism showing that dopamine denervation produces excessive 20–40 Hz synchronization of basal ganglia spiking with motor cortical LFP oscillations (Walters et al., 2007; Mallet et al., 2008; Galati et al., 2009; Brazhnik et al., 2012). The strong involvement of primary motor cortex in abnormal synchronization in humans has therapeutic implications, supporting the possibility of interrupting pathological synchronization at the level of the cortex (Drouot et al., 2004; Wu et al., 2007; Li et al., 2012), or of using a cortically based signal detector as a driver of basal ganglia deep brain stimulation in a closed loop paradigm (Rosin et al., 2011).
STN spike-M1 gamma synchronization
Our most novel finding is that STN spiking is also synchronized to M1 broadband gamma activity. Spike time averaged cortical gamma activity has a highly distinctive pattern, with waves of alternating high- and low-gamma amplitude, whose peak amplitude typically precedes the timing of STN spikes (median latency of −48 ms). These gamma waves have an oscillatory envelope in the theta, alpha, or beta ranges, demonstrating an interaction between the phase of low-frequency rhythms and the amplitude of broadband gamma.
This finding has important implications for the role of cross-frequency interactions in movement disorders. Previous studies of STN LFPs have called attention to the possibility of phase–amplitude interactions as a biomarker of PD (Lopez-Azcarate et al., 2010; Ozkurt et al., 2011). These studies showed coupling between the beta rhythm and a narrow bandwidth, very high-frequency activity (250–300 Hz). Subsequently, we used M1 LFP recordings to demonstrate that PD patients have excessive coupling of the phase of the beta rhythm to broadband gamma amplitude (de Hemptinne et al., 2013). The present study directly links this abnormal cortical phase amplitude coupling to basal ganglia discharge, and indicates that epochs of M1 phase amplitude coupling statistically predict the timing of STN spiking. Further, we show a symptom and disease-specific pattern of spike-synchronized gamma activity: tremor severity is inversely related to the magnitude of spike-synchronized gamma activity in PD (at any frequency) (Fig. 4F; Table 3), and beta phase modulation of spike-synchronized gamma activity occurred uniquely in the parkinsonian state, not in craniocervical dystonia patients. In the cortex, we did not observe any narrow bandwidth, very high-frequency peaks of gamma activity in the 250–300 Hz range that were distinct from the broadband gamma range, as has been observed in the STN LFPs (Foffani et al., 2003; Lopez-Azcarate et al., 2010; Ozkurt et al., 2011).
In classical models of PD, STN hyperactivity is thought to result from abnormal synaptic transmission through intrinsic basal ganglia pathways involving the striatum, primarily the substantia nigra compacta (SNc)-striatum-globus pallidus externa (GPe)-STN pathway (Bergman et al., 1990). Our findings suggest an additional mechanism for STN hyperactivity in PD. Because cortical broadband gamma activity is thought to reflect underlying action potential firing (Manning et al., 2009), we propose that M1 drives STN hyperactivity via the corticosubthalamic pathway (Nambu et al., 2002). The organization of M1 spiking activity into phase-modulated waves of high and low amplitude is ideally suited to drive STN activity, because strong phasic inputs to STN at beta envelope frequencies would drive STN spiking more efficiently than an asynchronous increase in corticosubthalamic activity (Beurrier et al., 1999; Plenz and Kital, 1999; Bevan et al., 2002), and at a faster rate than phasic inputs at lower envelope frequencies. The latency between M1 population spiking and STN spiking is longer than can be explained by fast monosynaptic neurotransmission. This may relate to the membrane dynamics of STN cells, in which repetitive slow depolarization is especially effective for triggering spikes in “burst mode,” thought to be characteristic of parkinsonism (Beurrier et al., 1999). A primary role for pathological activity in the corticosubthalamic pathway is also attractive in light of recent investigations in rodent models of PD, showing that modulation of the corticosubthalamic pathway is important for amelioration of motor deficits (Gradinaru et al., 2009; Li et al., 2012). Further, modeling work indicates that network oscillations in the basal ganglia-thalamocortical circuit are enhanced in the setting of excessive corticosubthalamic driving of intrinsic basal ganglia pathways (Holgado et al., 2010).
Limitations
All recordings studied here were performed in a state of alert rest. We did not study patients during a movement task due to the extreme difficulty of holding stable unit activity for a long period of time during patient movement. We did not study patients in the on medication state, so it is not clear how dopamine replacement therapy would affect the observed STN-cortex interactions (Williams et al., 2002). Due to the smaller number of dystonia subjects implanted, the number of units in dystonia patients with sufficient spike numbers for comparison with PD was relatively small. Because our PD subjects had longstanding disease and had received chronic levodopa, our findings may represent compensatory changes rather than primary abnormalities, or may reflect long-term effects of chronic levodopa treatment that do not “wash out” in 12 h off of levodopa (Picconi et al., 2003).
Conclusions
Using the novel approach of simultaneous recording of M1 LFPs and STN single units in awake humans undergoing DBS implantation surgery, we studied basal ganglia–cortex interactions on a fast time scale. In humans, we confirm the evidence from animal models that basal ganglia unit discharge in the parkinsonian state is strongly synchronized to cortical oscillatory activity. We further show that STN spiking is synchronized with cortical broadband gamma, and that the latter occurs in a phase-modulated pattern that begins before the occurrence of STN spikes. Abnormal cortical phase-amplitude coupling may be an important mediator of network oscillations in the basal ganglia-thalamocortical circuit.
Footnotes
This work was supported by the National Institutes of Health (Grant R01NS069779 to P.A.S.), and the Dystonia Medical Research Foundation. We thank Leslie Markun and Marta San Luciano for assistance with clinical data and statistical analyses, and Coralie de Hemptinne for critical review of the paper.
The authors declare no competing financial interests.
- Correspondence should be addressed to Dr. Philip Starr, Department of Neurological Surgery, University of California, 779 Moffitt, 505 Parnassus Avenue, San Francisco, CA 94143. starrp{at}neurosurg.ucsf.edu
