Abstract
Parkinson's disease (PD) is a heterogeneous disorder that leads to variable expression of several different motor symptoms. While changes in firing rate, pattern, and oscillation of basal ganglia neurons have been observed in PD patients and experimental animals, there is limited evidence linking them to specific motor symptoms. Here we examined this relationship using extracellular recordings of subthalamic nucleus neurons from 19 PD patients undergoing surgery for deep brain stimulation. For each patient, ≥10 single units and/or multi-units were recorded in the OFF medication state. We correlated the proportion of neurons displaying different activities with preoperative Unified Parkinson's Disease Rating Scale subscores (OFF medication). The mean spectral power at sub-beta frequencies and percentage of units oscillating at beta frequencies were positively correlated with the axial and limb rigidity scores, respectively. The percentage of units oscillating at gamma frequency was negatively correlated with the bradykinesia scores. The mean intraburst rate was positively correlated with both bradykinesia and axial scores, while the related ratio of interspike intervals below/above 10 ms was positively correlated with these symptoms and limb rigidity. None of the activity parameters correlated with tremor. The grand average of all the significantly correlated subthalamic nucleus activities accounted for >60% of the variance of the combined bradykinetic-rigid and axial scores. Our results demonstrate that the occurrence of alterations in the rate and pattern of basal ganglia neurons could partly underlie the variability in parkinsonian phenotype.
Introduction
Abnormal activity of basal ganglia neurons is thought to underlie the motor symptoms of Parkinson's disease (PD). The subthalamic nucleus (STN) has received particular attention in this regard, given the possibility of recording from this structure during the implantation of deep brain stimulation (DBS) electrodes.
It has been proposed that widespread changes in the firing patterns of basal ganglia neurons could be a key pathophysiological mechanism in PD. Both STN units and local field potentials (LFPs) exhibit oscillations in the sub-beta, beta, and gamma frequency ranges (Steigerwald et al., 2008; Eusebio and Brown, 2009; Schrock et al., 2009). The reduction in the power of beta oscillations following dopamine-replacement medication and DBS positively correlates with clinical improvement (Kühn et al., 2006, 2008, 2009; Ray et al., 2008). However, the proportion of STN neurons oscillating at beta frequency does not directly correlate with OFF clinical scores (Weinberger et al., 2006; Steigerwald et al., 2008; Zaidel et al., 2010). More recent studies have demonstrated that complex measures of beta oscillations in the LFP are weakly correlated with the severity of baseline symptoms (Chen et al., 2010; Pogosyan et al., 2010; Little et al., 2012), but the lack of a clear association with specific symptoms remains a confound to establishing them as causally pathophysiological (Stein and Bar-Gad, 2013). In PD patients and dopamine-depleted animals, STN units also display more bursting activity than in other neurological diseases and control animals, respectively (Bergman et al., 1994; Levy et al., 2000; Steigerwald et al., 2008; Tachibana et al., 2011). The intraburst rate in STN neurons is a highly accurate classifier of healthy and PD animals (Sanders et al., 2013), but there is no evidence that the extent of bursting is correlated with the clinical phenotype in PD patients.
Alternatively, the Albin–DeLong model of basal ganglia function suggests that the loss of dopaminergic input to the striatum results in a series of alterations in firing rate across cortico-basal ganglia-thalamic circuits that lead to a deficit in the initiation of movement (Albin et al., 1989; DeLong, 1990). Although there has been much dispute around some aspects of this model, it successfully predicts that the firing rate of STN units increases up to twofold in PD (Bergman et al., 1994; Mallet et al., 2008; Steigerwald et al., 2008; Schrock et al., 2009; Tachibana et al., 2011). However, STN firing rate does not correlate with motor symptom severity (Steigerwald et al., 2008).
Despite the abundance of evidence that changes in the activity of basal ganglia neurons are associated with the parkinsonian brain (Sanders et al., 2013), the relationship of different activities to specific PD symptoms is unclear. Here we examine several activity parameters of STN units recorded during the implantation of DBS electrodes in PD patients and retest the hypothesis that they can predict motor symptom severity.
Materials and Methods
Patient details and clinical scores.
This study was conducted in agreement with the Code of Ethics of the World Medical Association (declaration of Helsinki, 1967) and was approved by the local ethics committee. All patients gave their written informed consent to participate in this study. We studied 19 patients (7 female, 12 male; mean age ± SD: 65 ± 5 years) suffering from advanced idiopathic PD (average Hoehn and Yahr score ± SD: 3 ± 0.6; Hoehn and Yahr, 1967) with a disease duration of 13 ± 5 years (mean ± SD). The daily levodopa-equivalent dose was 1399 ± 566 mg (mean ± SD, conversion factors as described previously; Tomlinson et al., 2010). All patients met the selection criteria stated by the CAPSIT-PD committee (Defer et al., 1999) and underwent bilateral microelectrode-guided implantation of DBS electrodes in the STN. Preoperatively, all patients had significant improvement of the Unified Parkinson's Disease Rating Scale (UPDRS; Fahn et al., 1987) motor section (part III) score following intake of levodopa (mean symptom improvement ± SD: 51 ± 15%). Furthermore, all patients had an adequate global intellectual capacity [Mattis Dementia Rating Scale (Mattis, 1988), mean score ± SD: 139 ± 5], no structural alterations on MRI, and no concomitant severe medical comorbidities. Further clinical details are summarized in Table 1. UPDRS III scores were assessed 1 week (5 ± 3 d) before the DBS surgery, both OFF and ON medication (mean motor score in Dopa-OFF condition ± SD: 38 ± 11 vs Dopa-ON: 19 ± 8; p < 0.001, paired t test). Four subscores (OFF medication) were used for the present analysis: limb rigidity (the sum of UPDRS III sub-item 22 for arm and leg), bradykinesia (the sum of UPDRS III sub-items 23–26 and 31), axial score (the sum of UPDRS III sub-items 27–30), and tremor (the sum of UPDRS III sub-items 20 and 21).
Surgical procedures.
Dopamine agonist treatment was stopped >7 d before the operation, and surgery was performed after overnight withdrawal of anti-parkinsonian medication. Operations were performed under local anesthesia. Details concerning the surgical procedure were reported previously (Hamel et al., 2003). Briefly, an MRI-compatible Zamorano–Dujovny frame (Stryker Leibinger) was mounted on the patient's head and tightly secured with pins. Both gadolinium-enhanced volumetric T1 MRI and T2-weighted spin echo MRI sequences were acquired, and were fused with a computerized tomography scan using a commercially available algorithm (iPlan; Brainlab). Both commissures, the STN/nigra complex, and vessels were delineated with high resolution. After determining a reference line connecting the anterior and posterior commissure (AC–PC line), the STN was targeted 11–13 mm lateral to the midline, 1–3 mm inferior, and 1–3 mm posterior to the mid-commissural point on both sides. All trajectories used for microelectrode recording avoided blood vessels, sulci, and ventricles. A burr hole was fashioned anterior to the left and right coronal suture under local anesthesia, the micromanipulator was mounted on the stereotactic frame, and the appropriate target coordinates were adjusted. Ten to 15 minutes before the start of microelectrode recordings at thalamic level (typically 6–12 mm above the center of the STN as delineated on MRI), systemic analgesia with low-dose remifentanil (0.01–0.1 μg/kg/min) was completely stopped. No sedatives or anesthetics were administered.
Microelectrode recordings.
All intraoperative recordings (A.S., A.G., and C.K.E.M.) and subsequent signal analyses (A.S) were performed without prior knowledge of patients' presurgical symptom severity, which was assessed independently by movement disorder neurologists (S.Z. and C.B.). Microelectrode recordings were performed along five parallel tracks arranged in a concentric array (MicroGuide; Alpha Omega). In this configuration (originally introduced by A.-L. Benabid, “BenGun”), four outer platinum–iridium electrodes (mean impedance ± SD = 0.7 ± 0.2 MΩ at 1000 Hz; FHC) were separated by 2 mm from a central one, which aimed at the theoretical target. Signals were amplified (×20.000), bandpass-filtered (300–6000 Hz), and digitized (sampling rate: 24 kHz). Unit activities were identified as being from the STN using several previously described criteria (Sterio et al., 2002; Hutchison et al., 2004; Shamir et al., 2012). An electrode was identified as having entered the STN following a clear increase in the size of the background activity (Moran et al., 2006). STN units were distinguished by their tonic irregular, oscillatory, or bursty discharge pattern (Fig. 1), which was clearly different from slower bursting and single spiking of neurons in the overlying thalamus/zona incerta and from the high-frequency regular spiking of substantia nigra pars reticulata neurons found ventrally (Sterio et al., 2002).
Reconstructed stereotactic coordinates.
The stereotactic positions of the microelectrode tips for each unit recording were reconstructed using StereoTaxi, an in-house generated software. The position of the electrode tip relative to the mid-commissural point of the patient was computed as the vector a′ with a1′, a2′, and a3′ reflecting the lateral, frontal, and dorsal direction of the patient's coordinate system, respectively, as follows: Starting with a virtual trajectory being aligned to the dorsal direction axis of the patient, a reflects the vector pointing from the mid-commissural point of the patient toward the tip of the electrode of interest. In this scenario, a1 and a2 are the coordinates of this vector pointing from the central electrode of the BenGun configuration toward the electrode of interest, and a3 equals the distance between the electrode tip and its target position. Next, the BenGun configuration might be rotated along the longitudinal axis of the trajectory, which is reflected by the rotation matrix RBG. The virtual trajectory is then rotated according to the frontal and sagittal angles of the stereotactic frame and this second rotation is represented by the matrix RT. Finally, the shift between the mid-commissural point and the target position is added by vector t. These coordinates were normalized to the entrance of the STN (as identified physiologically) by subtracting the coordinates of the first position at which any microelectrode entered the nucleus from all other coordinates where an STN unit was recorded. We also performed a 2D normalization to the entrance of the STN whereby each position was given in relation to the entrance of the individual microelectrode (i.e., the position in 3D space is determined by the trajectory of the individual microelectrode) to the STN, as used in other studies (Weinberger et al., 2006).
Data processing.
All data were collected as part of our routine neuronavigation process, which aims to identify precisely the surgical target (STN). Patients were only included in this study if a minimum of 10 STN unit recordings, each lasting >45 s, were recorded from one or both hemispheres. The range of the number of units for each patient was 10–41 (mean = 21.2 ± 9.8; Table 1). This allowed us to calculate meaningful averages and SEMs within patients, rather than across the whole population of recorded units. It is important to note that recordings were only included if they were made during sustained periods where patients were awake and highly cooperative (vigilance was continuously monitored throughout the recording process). In total, 402 STN unit recordings were used in this study. In some recordings (118/402 units recorded in 15/19 patients), patients were engaged in a simple, brief movement task (flexion/extension of wrist or ankle) during short periods of the recording as part of the mapping procedure.
Spike detection was performed off-line using a voltage threshold method, which was set sufficiently high relative to the noise level to avoid false positives (only threshold values >4 SD of the background level were used). When possible, single units were then separated by manual cluster selection in 3D feature space on the basis of several waveform parameters including principal components, signal energy, and, importantly, the presence of a central trough in the autocorrelogram (Offline Sorter; Plexon). Over half of the unit activities were characterized as single-unit activities (SUAs; n = 287). Given the relatively high density of neurons in the STN, we chose to be conservative in spike sorting and in the classification of single units. Many recordings classified as multi-unit activities (MUAs) probably, therefore, comprised mostly one or two STN units and were clearly distinct from background spiking activity used in some previous studies (Moran et al., 2008; Zaidel et al., 2010). These multiple units were also present throughout the extracted/analyzed periods; they just could not be fully separated using our spike-sorting methods. In all cases, the initial 5–15 s of the unit recording was discarded as were any portions of injury discharge at the end of the recording.
Signal processing.
Spectral parameters of units were evaluated using fast Fourier transform as described previously (Halliday et al., 1995) using 0.5 Hz frequency bins. A Hanning window filter was used for all spectral analyses and spectra were estimated by averaging across these discrete sections (Halliday et al., 1995). Spectral analyses and related correlations were computed with a combined population of SUAs and MUAs and were calculated using the MATLAB toolbox NeuroSpec (www.neurospec.org). Significance was evaluated using surrogate data where the spike train was reconstructed by shuffling the interspike intervals (ISIs) across the record (Rivlin-Etzion et al., 2006a; Sharott et al., 2009). Power spectra were calculated on 1000 of these surrogate units and then used to construct 95% confidence limits for each frequency bin. To detect unambiguously significant peaks, we only considered a unit to be oscillating at a given frequency when the power in three contiguous bins (1.5 Hz total) exceeded the 95% confidence limit within that frequency range. Analysis of shuffled, thresholded spectra from 2 to 90 Hz (where zeroes represent spectral values lower than the confidence limit and ones represent significant values exceeding the confidence limit), showed that the probability of three contiguous bins exceeding the confidence limit by chance was ∼0.00002. The combination of the 95% confidence limit and three contiguous bin thresholds gave, therefore, a highly conservative estimate of when a unit was oscillating at a given frequency band and significant peaks were highly unlikely to represent noise. We calculated the percentage of units with a significant peak in the sub-beta (2–13 Hz), beta (14–31 Hz), and gamma (40–90 Hz) ranges in each patient. To calculate the spectral power, we first averaged the power spectra across all units recorded in a single hemisphere. As power was often concentrated in one region of the frequency band, we calculated the maximum power (of the mean spectrum) in a given band and summed it together with the values of the three frequency bins on either side (3.5 Hz total). We normalized this value by expressing it as the percentage of the total power in the 2–90 Hz range (of the mean spectrum). For patients where both hemispheres were recorded, we averaged the values for the two hemispheres. To calculate correlations, we used the logarithm of this value.
The firing rate of each neuron was calculated in the most stable part of the spike train. STN units are autonomously active (Do and Bean, 2003), and their tonic firing does not generally fall below 10 spikes/s in well isolated recordings (Magill et al., 2000; Mallet et al., 2008). The rate of each neuron was calculated using the longest part of the spike train where the firing frequency was >10 spikes/s in every 2 s bin for >30 s. The mean rate of single units was 32.9 spikes/s, in agreement with well isolated units in previous studies (Levy et al., 2000; Schrock et al., 2009; Weinberger et al., 2006; Moran et al., 2008). Bursts were detected using the Poisson surprise (PS) method, which has been used previously to detect bursts in the STN of parkinsonian primates (Legéndy and Salcman, 1985; Wichmann and Soares, 2006; Sanders et al., 2013). For a recorded spike train T, the probability of a contiguous subset of spikes t (defined in time) was quantified by constructing a Poisson-distributed spike train with the same rate as T and determining the rate of occurrence of t in the Poisson-distributed train. The negative logarithm of this probability is termed PS. Higher PS values indicate more improbable events in Poisson probability space. Sequences of spikes with a PS of ≥3 were considered to be bursts (Sanders et al., 2013). From this analysis, the proportion of spikes that occurred during bursts and the intraburst rate (mean spikes in burst/length of burst) was calculated for each unit.
Correlation and regression.
Measures of the magnitude (e.g., mean spectral power) or prevalence (e.g., the percentage of units with a significant peak in the power spectra) of activity parameters for each patient were plotted against the clinical subscores (limb rigidity, bradykinesia, axial symptoms, and tremor) OFF medication. When both variables were normally distributed (Shapiro–Wilk Test, p ≤ 0.05, SigmaPlot 12; Systat Software) and the variables appeared to be linearly related, a bivariate regression model was used to calculate the Pearson correlation coefficient (R, R2 values are shown) and associated p value (pr) to evaluate the strength, significance, and the best linear fit of the correlation. These criteria were met for all pairs of variables other than those including the tremor score, where the Spearman ρ coefficient was calculated (primarily due to the non-normal distribution of tremor scores). Where Pearson's correlation coefficient was used, values with a Cook's distance of > 1 were not included in the regression model (Cook, 1977). Values that met this criterion were only present in a small minority of correlations, and none of the significantly correlated variables had such outliers. In addition, we calculated a second p value using permutation statistics to test the hypothesis that the R value for the correlation between a given STN activity parameter and clinical score could result from a random distribution of the physiological data, rather than being linked to the specific distribution across these patients (and associated clinical scores). To achieve this, we created null-hypothesis sample datasets by shuffling the values for the STN activity parameters across the whole-unit sample (402 all units/271 SUAs). The STN activity parameter variables for each patient were then calculated using the same number of these shuffled units that was recorded for each patient in the real data (i.e., Table 1, last column) and the clinical scores were not changed. Thus, the null-hypothesis samples differed only in the distribution of the real values for each unit across the patient sample. The R value was calculated for the correlation between 10,000 of the null-hypothesis samples to give a null R value distribution. The p value for this permutation statistic (pp) was the probability of the real R value (from regression model of the real data) occurring in the null R value distribution. Thus, a pp of 0.05 indicates that there is a 5% chance that the R value would occur in the distribution of R values calculated using the null samples.
A total of 40 pairs of STN activity parameter and clinical score were tested (Table 2) using these methods and were corrected for multiple comparison using the false discovery rate (FDR) method (Benjamini and Hochberg, 1995). Correlations that had pr and pp values < 0.05 only before FDR corrections were designated as “weakly correlated.” Correlations that had significant pr and pp values following FDR correction were designated as “significantly correlated.” We describe the weakly correlated values to provide more detail about the selectivity of correlations between given STN activity parameters and different symptoms. Multiple-comparison correction was performed only on the variables used to test the original hypotheses. Post hoc correlations to refine details or adjustments of these core observations (e.g., laterality, spatial bias, effect of movement periods) were performed only on variable pairs that were significant in the main analysis and were considered significant at the p < 0.05 level. In some cases, partial correlations were used to determine the dependence of a correlation on a third variable and were considered significant at p < 0.05.
Correlations were also calculated using indices comprised of multiple STN activity parameters. These analyses were performed only in the 16 patients that had ≥9 SUAs so that oscillation and ISI ratio parameters could be combined. The 16 values for each STN activity parameter were normalized using a z-score ([values for each patient − mean all patients]/SD all patients). These normalized scores (for single STN-activity parameters) could then be averaged (giving one set of 16 values where each STN-activity parameter was weighted equally) and regression models with the UPDRS subscores calculated. This process, therefore, involved no optimization to find the best-fit model for a given symptom, but allowed us to test whether using a simple combination of two or more activity parameters could improve correlations with individual or multiple symptoms, as compared with those parameters alone.
All statistical and signal analyses were performed using MATLAB (MathWorks) and associated toolboxes. All values are given as means ± SEM unless noted otherwise.
Results
The spontaneous activity of 402 STN units was recorded from 19 PD patients during the mapping of the subthalamic area. STN units displayed clear variation in their firing rate and their propensity to burst and oscillate, within and between patients (Fig. 1).
Linear regression analysis of STN activity parameters and motor symptom severity
We examined the correlation between four parkinsonian motor symptoms (limb rigidity, bradykinesia, axial symptoms, and tremor), as measured by the preoperative UPDRS OFF scores, and 10 STN activity parameters. For three frequency bands (sub-beta, beta, and gamma) we made separate calculations for oscillatory power and the prevalence of oscillating units (regardless of the strength of the oscillations), resulting in six activity parameters that measured neuronal oscillations. The remaining four activity parameters measured different aspects of firing rate and bursting. This gave a total of 40 comparisons (Table 2). For each pair of STN activity parameter and symptom, we calculated a p value for both the regression model and for the distribution of the underlying data. After performing an FDR correction of multiple comparisons on both sets of p values, eight STN activity parameter/symptom pairs were found to be significantly correlated. Using the p values for Spearman correlations with the non-normally distributed tremor scores did not change these results. In the following section, we demonstrate these significant correlations and describe the relevant STN activity parameters in detail.
STN activity parameters correlate with the severity of specific parkinsonian motor symptoms
STN units were examined for oscillatory activity from 2 to 90 Hz. Clear peaks in the power spectra exceeding the 95% confidence intervals could be observed for many individual units (Fig. 2A). As reported previously (Weinberger et al., 2006), many neurons displayed prominent peaks in the sub-beta and beta frequency ranges (Fig. 2A). Sharper peaks were also observed at gamma frequencies in individual units (Fig. 2B). There were clear peaks in the percentage of significantly oscillating units across patients at sub-beta and beta frequencies (Fig. 2C).
We examined whether the mean spectral power across all units and proportion of units oscillating at sub-beta (2–13 Hz), beta (14–31 Hz), and gamma (40–90 Hz) frequencies correlated with the preoperative UPDRS subscores for all patients (n = 19). The mean spectral power in the sub-beta band was weakly correlated with the bradykinesia score (Table 2) and significantly positively correlated with the axial score (Fig. 2D). Dividing the sub-beta band into theta and alpha bands did not lead to significant correlations. The percentage of units oscillating at beta frequency was weakly positively correlated with the bradykinesia score (Table 2) and significantly positively correlated with the limb rigidity score (Fig. 2E). The proportion of units oscillating at gamma frequency was strongly negatively correlated with the bradykinesia score (Fig. 2F). There were no significant correlations with the tremor subscores. In summary, the oscillatory firing of STN neurons in different frequency bands was selectively associated with different motor symptoms.
Next we examined the hypotheses that an increase in firing rate and/or bursting of STN units could underlie motor symptoms in PD patients. Each STN unit was classified as exhibiting SUA or MUA based on several criteria (see Materials and Methods). It is important to note that we only included SUAs when correlating firing rate and burst parameters with clinical scores, as these parameters become difficult to interpret when used on MUAs. Thus we removed three patients who had eight or less single units, which was substantially lower than the minimum number of units used to evaluate correlations with oscillation parameters (i.e., 10). The mean (32.9 spikes/s) firing rate of these SUAs (Fig. 4A) was in good agreement with well isolated units in parkinsonian patients (Levy et al., 2000; Moran et al., 2006; Weinberger et al., 2006; Schrock et al., 2009). Using this unit/patient sample (n = 16), we found a weak positive correlation between the mean firing rate and the axial scores (Table 2). Using values for mean firing rate calculated with no control for stability (see Materials and Methods) did not lead to any significant results.
Using the PS method, we calculated percentage of spikes in bursts for each unit and averaged this for each patient. This parameter was not significantly correlated with any of the clinical scores. Based on a recent study of features that distinguish STN neurons from healthy and parkinsonian primates (Sanders et al., 2013), we then calculated the correlations with the mean intraburst rate (number or spikes/burst length). This measure was weakly positively correlated with the limb rigidity score (Table 2) and significantly positively correlated with the bradykinesia (Fig. 3A) and axial (Fig. 3B) scores. Neither of the burst parameters was significantly correlated with the tremor score.
PS-based burst analysis has several parameters that must be arbitrarily selected (e.g., PS threshold). Moreover, STN neurons can display high instantaneous firing episodes both inside (Fig. 4Ai–Aiii) and outside (Fig. 4Aiv) of bursts. Thus, we examined whether a more simple measure of high instantaneous firing rate, the ratio of ISIs below/above 10 ms (ISI ratio), was also correlated with the clinical scores. Ten milliseconds was chosen because repetitive firing with this ISI would be outside our gamma frequency range for the spectral analysis. Furthermore, this ISI ratio has previously been used in physiology studies of the human basal ganglia (Favre et al., 1999). Unsurprisingly, this ISI ratio was significantly correlated to both mean firing rate and intraburst rate at both the neuronal and patient levels (Fig. 4B–D). However, the ISI ratio was significantly correlated with the limb rigidity, bradykinesia, and axial scores, with higher R2 values than the intraburst rate (Fig. 4F–H). The correlation between the ISI ratio and the sum of these scores was even stronger, accounting for >50% of their variance between patients. Together, these results suggest that the frequency of high instantaneous firing episodes is strongly and nonselectively associated with motor symptom severity. Measures of mean firing rate and intraburst rate may give diluted information about this underlying parameter. Thus, we used the ISI ratio to represent these high firing rate episodes for the remainder of the study.
For the STN activity parameter/clinical score pairs described above, similar p and R2 values were also observed using nonparametric Spearman correlation coefficients.
Bradykinesia and limb rigidity are often lateralized in PD patients. Thus, we re-examined the significant correlations involving these symptoms using units from single hemispheres and the contralateral and ipsilateral hemibody scores for those hemispheres. This analysis was conducted in the 25 hemispheres where ≥10 units were recorded. There were significant correlations (p < 0.05) for both ipsilateral and contralateral hemispheres and, therefore, no significant differences when considering hemispheres separately. This lack of lateralization may have been partly due to a strong positive correlation between the ipsilateral and contralateral bradykinesia and rigidity scores (R2 = 0.44, p = 0.0002), which would lead to the hemibody scores in each patient being predictive of each other. Given this apparent lack of lateralization, analysis of hemibody scores was not considered further. We also calculated a lateralized tremor score based on the occurrence of tremor during the operation (based on intraoperative EMG recordings and clinical observation: the percentage of units recorded with a tremor period in either contralateral limb). This measure was not significantly correlated with any of the STN activity parameters.
Movement episodes could not account for correlations between STN activity parameters and motor symptoms
During the recording of some units, patients were asked to make brief (5–20 s), sustained flexions and/or extensions of the ankle or wrist. This was necessary to map the motor/non-motor portions of the STN, but due to time constraints it was not always possible to make separate rest and movement recordings for each unit. Although the movement period usually represented only a small proportion of the recording, it could be a confounding variable in the symptom/STN activity correlations. The percentage of STN unit recordings in each patient containing a movement period, however, was not significantly correlated with any clinical scores (Pearson correlation, pr < 0.05). In addition, keeping the same constraints on the minimum sample (≥9 units) included for each patient (n = 14), the majority of the strongest correlations were still significant at the p < 0.05 level following the removal of all units that were recorded in periods where there was cued movement (axial/sub-beta power: R2 = 0.42, pr = 0.0126; limb rigidity/beta units: R2 = 0.4, pr = 0.015; bradykinesia/gamma units: R2 = 0.35, pr = 0.015). In the case of the ISI ratio, removing all single-unit periods with movement episodes halved the data for both total units and patients (187 units, 8 patients) and the correlations were no longer significant. However, correlations between the ISI ratio and rigidity (R2 = 0.47, p = 0.0048), bradykinesia (R2 = 0.37, p = 0.016), and axial symptoms (R2 = 0.3, p = 0.033) were still significant at the p < 0.05 level following partial correlation of the percentage of units recorded with a movement period. Despite substantial reductions in the amount of data used, the key findings could, therefore, all be observed following the removal of the influence of movement periods.
Bias in recording position could not account for correlations between STN activity parameters and motor symptoms
Previous studies suggest that beta oscillating units are more commonly found at the dorsolateral aspect of the STN (Weinberger et al., 2006), while gamma activities are stronger in more ventral areas (Zaidel et al., 2010). If these distributions were present in our data, a bias toward recording in these areas in patients of a given phenotype could lead to spurious correlations. To this end, we examined the location of the STN units in two ways: the 2D depth normalized to the entrance on each individual electrode (Weinberger et al., 2006), and the 3D coordinates normalized to the first electrode to enter the STN (Fig. 5). Using the normalized 3D coordinates, units oscillating in the beta band appeared to cluster toward the dorsolateral corner of the STN (Fig. 5A), but no significant difference was found in any plane (data not shown). When depth was examined on the 2D axis for each electrode, however, beta oscillating units were approximately twice as frequent between 0.5 and 1.5 mm from the STN entrance as at other depths (Fig. 5B) and were significantly closer to the entrance than nonoscillating units (Mann–Whitney, p = 0.0016; Fig. 5C). As the percentage of beta units was most strongly correlated with rigidity, we examined whether a bias in recording position could influence this correlation. The mean recording position was not significantly correlated with the rigidity score, suggesting we did not record more in areas with more beta units in patients with higher rigidity scores (Fig. 5D). In contrast, gamma oscillating units appeared to be more common in the medial STN (Fig. 5E,F) and were significantly further from the STN entrance along the normalized x-axis (mediolateral, Mann–Whitney, p = 0.0063; Fig. 5G). The percentage of gamma oscillating units was most strongly correlated with bradykinesia, but there was no correlation between the mean recording position on the x-axis and this score (Fig. 5H). Sub-beta power was significantly further from the STN entrance along the normalized dorsoventral z-axis (Pearson correlation; Fig. 5I). The sub-beta power was strongly correlated with axial score, but there was no correlation between the mean recording position on the z-axis and this score (Fig. 5J). No significant correlations were found between high-frequency firing and any axis. There were also no significant correlations between clinical scores and the percentage of units categorized as being in a specific part of the nucleus (e.g., dorsal vs ventral; data not shown). In summary, there were biases in the prominence of neuronal oscillations in different STN territories, but this bias was unlikely to effect the correlations between STN unit parameters and symptom severity.
Selectivity of STN activity parameters in predicting clinical subscores
Previous studies have demonstrated relationships between oscillations at different frequencies in STN LFPs recorded from PD patients (Fogelson et al., 2005; Marceglia et al., 2006). However, there was no significant correlation between the percentage of units oscillating at beta frequency and the percentage of units oscillating at gamma frequencies across patients. Thus, in contrast to studies of LFP power (Fogelson et al., 2005), the number of units oscillating at beta frequency was not related to the number of units oscillating at gamma frequency. At the level of units, only 5% of units with beta peaks had gamma peaks and 12% of units with gamma peaks had beta peaks. The sub-beta power was also not significantly correlated with either of these parameters. Generally, this concurs with each type of oscillation having a different spatial bias and being associated with different symptoms. The oscillatory parameters, therefore, displayed little redundancy, partly explaining why they showed selectivity in their correlations with clinical scores. In contrast, both the percentage of beta oscillating units (R2 = 0.39, p = 0.01) and the mean sub-beta power (R2 = 0.37, p = 0.012) were positively correlated with the ISI ratio. The percentage of gamma oscillating units was not correlated with the ISI ratio. Sub-beta power, beta oscillating units, and the ISI ratio were therefore associated to some degree, while gamma oscillating units were independent of all other parameters.
Having investigated individual STN activities, defined the correlations between them, and excluded possible confounds from spatial bias and movement, we hypothesized that combinations of these parameters could lead to stronger and more specific correlations with motor symptoms than when compared alone. Thus, we calculated averages of the normalized percentages of the neurons displaying combinations of STN activity parameters and correlated them with UPDRS subscores (Fig. 6). Combining the beta-oscillating units and ISI ratio did not markedly improve on the predictive value of the limb rigidity for either of these variables alone (Fig. 6A). Similarly, combining the sub-beta power and ISI ratio led to only a modest improvement in predicting the axial score (Fig. 6B). In contrast, combining gamma oscillating units and the ISI ratio accounted for ∼20% more of the variance in the bradykinesia scores (Fig. 6C). This may therefore reflect the independence of the gamma oscillating units and ISI ratio, which could be correlated to different aspects of the bradykinesia scores. Combining all variables that were significantly correlated with limb rigidity, bradykinesia and axial scores accounted for nearly 65% of the variance of those scores combined (Fig. 6D). In summary, the firing pattern of STN units contains considerable information about motor symptom severity and disease phenotype.
Discussion
In patients undergoing implantation of DBS electrodes, we re-addressed the question of whether parameters of STN firing rate and pattern correlate directly with the severity of different motor symptoms in PD. In contrast to several previous studies, we were able to demonstrate significant correlations with the severity of motor symptoms.
STN activity parameters are correlated with the baseline severity of motor symptoms in PD
Dopamine depletion undoubtedly leads to changes in the firing rate and pattern of neurons in the basal ganglia (Hammond et al., 2007; Rosin et al., 2007). It has been recently demonstrated that these features are effective for classifying whether basal ganglia neurons are from a healthy or parkinsonian brain (Sanders et al., 2013). Despite this, few studies have demonstrated a direct correlation between basal ganglia unit activity and the severity of motor symptoms. Hitherto, correlations have been made either with the dopamine response (Kühn et al., 2005; Weinberger et al., 2006; Ray et al., 2008) or using complex measures of the LFP (Pogosyan et al., 2010; Yoshida et al., 2010).
Why did we find significant results where others have not? Uncovering correlations between neurological and physiological parameters depends on the signal-to-noise ratio of the clinical scores and electrophysiological data. We only included patients from whom we recorded 10 or more STN units, which allowed a robust calculation of each activity parameter. From the clinical standpoint, dopaminergic agonist treatment was withdrawn at least 1 week before the operation. Recordings were, therefore, made in the full OFF state and thus likely mirrored the clinical phenotype to a reasonable degree, despite possible sources of noise, such as microlesions. The large amount of data included for each patient may also have allowed us to detect significant correlations through this interference. To guard against false positive results, correlations were verified using permutation statistics and had to pass an FDR correction. Similar studies have not found significant results with OFF UPDRS scores, despite not applying such rigorous controls (Weinberger et al., 2006; Steigerwald et al., 2008). These findings are, therefore, likely the result high of signal-to-noise ratio of our dataset, rather than the result of performing extensive signal and statistical analysis.
Do STN unit activities cause parkinsonian symptoms?
The methods used here cannot be used to infer that these STN activities cause motor symptoms, particularly given that electrophysiological recordings and clinical evaluation were not coincident. Instead, it seems likely there is a strong association between the amount of underlying cellular pathology, the resulting level of abnormal basal ganglia activity, and the subsequent motor deficits (Jenkinson and Brown, 2011). Validation of the motor UPDRS suggests it is robust over the timescale of weeks (Siderowf et al., 2002). Abnormal basal ganglia activity is also stable over these time periods, as evidenced by the consistency of frequency profiles of STN-LFPs recorded weeks apart (Bronte-Stewart et al., 2009). Adequate electrophysiological recording should therefore allow predictions about the preoperative symptoms, even when not coincident with their assessment. This principle is perhaps best illustrated by axial symptoms, which may not be evident and/or measureable during surgery. However, if many units are recorded over a relatively long time period, the underlying axial symptom-associated STN activities should be more common in patients with severe symptoms than those without.
The wider question of causality is complicated by several other factors. First, symptom-related electrophysiological activities can be secondary to those in the true disease locus, rather than driving behavior themselves. For example, it is common to encounter STN units with firing that is strongly synchronized to the tremor (Rodriguez et al., 1998), but these are likely to be the result of afferent inputs (Rivlin-Etzion et al., 2006b). Second, symptom-associated activities could be the result of compensatory changes that increase in parallel with the primary deficit (Fan et al., 2012), which could explain the lag between dopamine depletion and the increase of oscillations (Leblois et al., 2007). However, for a given neuronal activity to cause a disease symptom, it could be considered a prerequisite that it predicts the severity of that symptom. The lack of such correlations has brought into question the role of STN beta oscillations as a pathophysiological mechanism in PD (Brown, 2007; Stein and Bar-Gad, 2013). Thus, while these correlations do not impart causality, they provide a stronger foundation for the candidacy of these STN activities as mechanistic factors.
Neuronal oscillations in the STN are selectively associated with motor symptoms
A major aim of this study was to investigate the selectivity of associations between STN activities and individual symptoms. The correlation among UPDRS scores complicates this issue (Martínez-Martín et al., 1994), but their use is necessary given that continuous assessment of phenotype during neurosurgery is impossible. Despite the redundancy in UPDRS items measuring rigidity, bradykinesia, and axial symptoms, we found some selectivity in their association with oscillatory activities. We provide the first demonstration of a positive correlation between STN sub-beta frequency oscillations and axial symptoms. Interestingly, the power of LFP oscillations in a similar range (6–10 Hz) in the pedunculopontine nucleus is also correlated with gait disturbance (Thevathasan et al., 2012). The finding that the proportion of beta oscillating units correlates positively with limb rigidity adds to the considerable evidence that these activities may have a mechanistic role (Jenkinson and Brown, 2011). In contrast, the number of units oscillating at gamma frequencies was negatively correlated with bradykinesia. This finding is in line with the positive correlation between gamma power and coherence with symptom improvement following dopamine replacement (Litvak et al., 2012). The percentage of gamma oscillating units was the most independent activity parameter and therefore the best candidate for being selectively associated with an individual symptom: bradykinesia. Together, these results are consistent with the hypothesis that beta oscillations in the motor system reinforce the status quo (Gilbertson et al., 2005; Engel and Fries, 2010), while gamma oscillations facilitate a readiness to move (Schoffelen et al., 2005).
In contrast, the intraburst rate and frequency of high-firing periods showed little or no selectivity among the nontremor symptoms, but accounted for >50% of the variance of these symptoms together. This supports the finding that the STN intraburst frequency was most predictive of the parkinsonian state in MPTP-treated primates (Sanders et al., 2013). Consequently, the increase in STN discharge in PD could be secondary to the increases in high instantaneous firing episodes. This is supported by the finding that the extraburst firing rate is not predictive of the parkinsonian state (Sanders et al., 2013). High instantaneous firing episodes in STN need not necessarily cause network disruption to predict disease severity. Indeed, therapeutically effective high-frequency DBS can entrain basal ganglia neurons (Bar-Gad et al., 2004; Meissner et al., 2005), which would generally reduce ISIs. Rather, an increased incidence of high firing rate episodes could be the result of increased entrainment to specific phases of network oscillations (Moran et al., 2008), which would reduce the temporal window in which spikes occur. In this case, such activities could provide a proxy for abnormal network synchronization overall, which could explain their positive correlation with several motor symptoms. As previously suggested, rather than modulating firing rate, the therapeutic effect of DBS could result from reshaping the temporal structure of neuronal activity (Bar-Gad et al., 2004; Xu et al., 2008), dissociating it from and thus disrupting pathological network oscillations (Meissner et al., 2005; Kühn et al., 2008).
Despite its relative independence (Martínez-Martín et al., 1994), none of the activity parameters analyzed here was correlated with the tremor score. In concordance with previous reports (Rodriguez et al., 1998; Rodriguez-Oroz et al., 2001), we observed STN units that fired in time with limb tremor (data not shown), but they did not predict the magnitude of the tremor score. This latter interpretation would be broadly consistent with the proposal that tremor severity has a predominantly cerebellar-based pathology (Stochl et al., 2008; Zaidel et al., 2009). However, we cannot rule out that more accurate characterization of tremor and/or different neurophysiological parameters would yield significant correlations, as shown in a recent study demonstrating a weak correlation between the spike-LFP modulation index and the contralateral tremor score (Shimamoto et al., 2013).
Implications for treatment
The combination of STN activity parameters studied here accounted for 64% of the variance on nontremor motor symptoms between patients, making STN activity nearly twice as predictive of the severity of these symptoms than the level of dopamine (Pirker, 2003). The discovery of accurate neuronal indicators of symptom severity has assumed greater clinical significance with the advent of closed-loop DBS, which has the potential to improve DBS efficacy (Tass et al., 2003; Rosin et al., 2011; Little et al., 2013). Long-term recording of single units has proven to be difficult in brain–computer interface applications; the activity parameters here may be most useful in helping to define the best “proxy” population activities for detecting disease-related activity “on the fly.” Given that symptoms change over time, with freezing and gait disruption becoming more common in later stages, using symptom-selective feedback signals has the potential to enhance the efficacy of such applications.
Footnotes
This work was supported by grants from the European Union (MRTN-CT-2005-019247 to A.K.E, FP7-ICT-270212 to A.K.E. and A.G., and Marie Curie European Re-integration Grant SNAP-PD to A.S.; Medical Research Council, UK (award U138197109, A.S.); and the Deutsche Forschungsgemeinschaft (SFB 936 to A.M., C.G., and A.K.E.). We are grateful to Peter Brown for valuable scientific discussions.
The authors declare no competing financial interests.
- Correspondence should be addressed to either of the following: Dr. Andrew Sharott, Medical Research Council Anatomical Neuropharmacology Unit, University of Oxford, Oxford OX1 3TH, United Kingdom, andrew.sharott{at}pharm.ox.ac.uk; or Dr. Christian Moll, Department of Neurophysiology and Pathophysiology, University Medical Center Hamburg-Eppendorf, 20246 Hamburg, Germany, c.moll{at}uke.uni-hamburg.de