Abstract
Tremor, a common and often primary symptom of Parkinson's disease, has been modeled with distinct onset and maintenance dynamics. To identify the neurophysiologic correlates of each state, we acquired intraoperative cortical and subthalamic nucleus recordings from 10 patients (9 male, 1 female) performing a naturalistic visual–motor task. From this task, we isolated short epochs of tremor onset and sustained tremor. Comparing these epochs, we found that the subthalamic nucleus was central to tremor onset, as it drove both motor cortical activity and tremor output. Once tremor became sustained, control of tremor shifted to cortex. At the same time, changes in directed functional connectivity across sensorimotor cortex further distinguished the sustained tremor state.
SIGNIFICANCE STATEMENT Tremor is a common symptom of Parkinson's disease (PD). While tremor pathophysiology is thought to involve both basal ganglia and cerebello–thalamic–cortical circuits, it is unknown how these structures functionally interact to produce tremor. In this article, we analyzed intracranial recordings from the subthalamic nucleus and sensorimotor cortex in patients with PD undergoing deep brain stimulation surgery. Using an intraoperative task, we examined tremor in two separate dynamic contexts: when tremor first emerged, and when tremor was sustained. We believe that these findings reconcile several models of Parkinson's tremor, while describing the short-timescale dynamics of subcortical–cortical interactions during tremor for the first time. These findings may describe a framework for developing proactive and responsive neurostimulation models for specifically treating tremor.
- deep brain stimulation
- directed functional connectivity
- electrocorticography
- Parkinson's disease
- subthalamic nucleus
- tremor
Introduction
Tremor, a cardinal symptom of Parkinson's disease (PD), typically manifests as a 4–6 Hz oscillatory movement of the distal limbs during rest or sustained posture (Lance et al., 1963). While often the presenting motor symptom of PD, tremor (and its response to dopamine replacement therapy) is highly variable across patients (Koller, 1984, 1986; Zach et al., 2015; Dirkx et al., 2017, 2019). PD tremor neurophysiology has been described by the “dimmer switch” model where an “on-off” mechanism is separable from a magnitude controller (Helmich et al., 2012). Specifically, functional magnetic resonance imaging (fMRI) BOLD activity from basal ganglia nuclei such as the globus pallidus pars interna (GPi) correlates with the presence or absence of tremor, whereas immediate tremor amplitude better correlates with BOLD signal from structures in cerebello–thalamo–cortical circuits such as motor cortex (Helmich et al., 2011; Helmich, 2018). The GPi and the monosynaptically connected subthalamic nucleus (STN; Albin et al., 1989) are common therapeutic targets for deep brain stimulation (DBS). Indeed, DBS in each nucleus is equally effective in reducing tremor (Wong et al., 2020). However, the precise role of the STN and its interactions with cortex in these tremor dynamics is unknown.
Low-frequency (4–8 Hz) oscillatory bursting has been observed in both the STN and GPi MPTP primate models of PD (Bergman et al., 1994; Raz et al., 2000). This bursting, although present in the absence of tremor, becomes highly synchronized with tremor once it emerges. STN recordings from patients with PD have similarly revealed theta/tremor frequency (3–8 Hz) activity that is coherent with electromyography (EMG) recordings of tremulous limbs (Levy et al., 2000; Reck et al., 2009, 2010). Accordingly, STN tremor frequency oscillations (along with higher-frequency oscillations) have been used to predict clinical measures of tremor (Hirschmann et al., 2016; Telkes et al., 2018; Asch et al., 2020). Further, studies applying STN DBS at tremor frequencies entrained tremor to the phase of the stimulation, consistent with a direct modulatory role of STN on tremor (Cagnan et al., 2014).
At the same time, tremor reorganizes cortical activity. Magnetoencephalography studies of patients with PD identified a broad cortical tremor network comprising “intrinsic” (ventrolateral anterior thalamus, premotor and motor cortex) and “extrinsic” (cerebellum, ventrolateral intermedius (VIM), somatosensory cortex) loops hypothesized to initialize and stabilize tremor, respectively (Volkmann et al., 1996; Timmermann et al., 2003). This corticocortical synchronization at single and double tremor frequencies extends to STN local field potential (LFP) and EMG recordings as well (Hirschmann et al., 2013). Meanwhile, intraoperative studies combining electrocorticography (ECoG) and STN LFP recordings found decreases in alpha (8–13 Hz) and beta (13–30 Hz) coherence during tremor (Qasim et al., 2016). Despite this broad corticocortical synchronization at tremor frequencies, it remains unclear whether these neurophysiological changes are specific to tremor onset or maintenance. In addition, although STN and sensorimotor cortex become coherent during tremor, the manner in which tremor-related activity is coordinated across structures and how different networks of activity may reflect the different stages of tremor production and maintenance are unknown.
Thus, to understand whether there are indeed distinct neurophysiological mechanisms of tremor initiation and maintenance, and to better understand what neurophysiological interactions characterize these states, we recorded local field potential activity from the STN along with ECoG from sensorimotor cortices while subjects with PD engaged in a task that elicited initiation and persistence of tremor. Specifically, we tested whether the STN (like the GPi) drove tremor specifically during onset, while cortical structures drove sustained tremor.
Materials and Methods
Participants.
All patients undergoing routine, awake placement of deep brain stimulating electrodes for intractable, idiopathic PD between November 2015 and September 2017 were invited to participate in this study. Patients with PD were selected and offered the surgery by a multidisciplinary team based solely on clinical criteria, and the choice of the target (STN vs GPi) was made according to each patient's particular circumstance (disease manifestations, cognitive status, and goals; Akbar and Asaad, 2017). In this report, we focused on 10 patients [9 male (M), 1 female (F)] undergoing STN DBS with intraoperative ECoG recordings. Patients were off all anti-parkinsonian medications for at least 12 h in advance of the surgical procedure (Unified Parkinson's Disease Rating Scale (UPDRS) part III score, 48.2 ± 15.6). Four patients were considered tremor dominant, and six patients had average tremor UPDRS III scores >2 in their right hand (Jankovic et al., 1990). Approximately age-matched control subjects [3 M, 11 F (often patients' partners)] also participated in this study (n = 14 subjects); patients were 55.6–78.5 years of age (mean age, 65.2 ± 7.4 years), and control subjects were 48.3–79.2 years of age (mean age, 62.4 ± 10.0 years) at the time of testing (Mann–Whitney U test comparing ages, p > 0.05). Control subjects were required simply to be free of any diagnosed or suspected movement disorder and to have no physical limitation preventing them from seeing the display or manipulating the joystick. There was a strong male bias in the patient population (9 M, 1 F) and a female preponderance in the control population (3 M, 11 F), reflecting weaker overall biases in the prevalence of PD and the clinical use of DBS therapy (Accolla et al., 2007; Hariz et al., 2011; Rumalla et al., 2018). All subjects were right handed. Patients and other subjects agreeing to participate in this study signed informed consent forms, and experimental procedures were undertaken in accordance with an approved Rhode Island Hospital human research protocol (Lifespan Internal Review Board Protocol #263157) and the Declaration of Helsinki.
Surgical procedure.
Microelectrode recordings (MERs) from the region of the STN of awake patients are routinely obtained to map the target area and guide DBS electrode implantation. A single dose of short-acting sedative medication (typically, propofol) was administered before the start of each procedure, at least 60–90 min before MERs. The initial trajectory was determined on high-resolution (typically, 3 T) magnetic resonance (MR) images coregistered with computed tomography (CT) images demonstrating previously implanted skull-anchor fiducial markers (version 3.0; FHC). Localization of the target relied on a combination of direct and indirect targeting, using the visualized STN as well as standard stereotactic coordinates relative to the anterior and posterior commissures. Appropriate trajectories to the target were then selected to avoid critical structures and to maximize the length of intersection with the STN. A 3-D printed stereotactic platform (STarFix Micro-Targeting System, FHC) was then created such that it could be affixed to these anchors, providing a precise trajectory to each target (Konrad et al., 2011). Microdrives were attached to the platform and then loaded with microelectrodes. Recordings were typically conducted along the anterior, center, and posterior trajectories (with respect to the initial MRI-determined trajectory) separated by 2 mm, corresponding to the axis of highest anatomic uncertainty based on the limited visualization of the STN on MRI. Bilateral ECoG strips were placed posteriorly along sensorimotor cortices through the same burr hole used for MER insertion for temporary recordings. MER began ∼10–12 mm above the MRI-estimated target, which was chosen to lie near the inferior margin of the STN, approximately two-thirds of the distance laterally from its medial border. The STN was identified electrophysiologically as a hyperactive region typically first encountered ∼3–6 mm above estimated target (Gross et al., 2006). At variable intervals, when at least one electrode was judged to be within the STN, electrode movement was paused to assess neural activity and determine somatotopic correspondence, as per routine clinical practice. At these times, if patients were willing and able, additional recordings were obtained in conjunction with patient performance of the visual–motor task.
Neurophysiological signals and analysis.
Microelectrode signals were recorded using “NeuroProbe” tungsten electrodes (Alpha Omega). ECoG signals were acquired using 8-contact subdural strips with 10 mm contact-to-contact spacing (Ad-Tech Medical). All signals were acquired at 22–44 kHz and synchronized using Neuro Omega data acquisition systems (Alpha Omega). Microelectrode impedances were typically 400–700 kΩ, while ECoG contact impedances were typically 10–30 kΩ. Patients performed up to four sessions of the task, with microelectrodes positioned at different depths for each session. As microelectrodes were not independently positionable, some signals may have necessarily been acquired outside of the STN. All recorded signals were nevertheless considered and analyzed.
Neural data were analyzed using the “numpy/scipy” Python 3 environment (Harris et al., 2020; Virtanen et al., 2020; https://numpy.org/, https://www.scipy.org/). Offline, ECoG contacts were rereferenced to a common median reference within a strip (Liu et al., 2015). All resulting signals were bandpass filtered between 2 and 600 Hz, and notch filtered at 60 Hz and its harmonics. Time series were z-scored and artifacts above 4 SDs were removed. These resulting time series were then downsampled to 1 kHz. Time series were bandpass filtered using a Morlet wavelet convolution (wave number 7) at 1 Hz intervals, covering 3–400 Hz. The instantaneous power and phase at each frequency was then acquired by the Hilbert transform. To analyze broad frequency bands, we grouped frequencies as follows: theta, 3–8 Hz; alpha, 8–12 Hz; low beta, 12–20 Hz; high beta, 20–30 Hz; low gamma, 30–60 Hz; mid-gamma, 60–100 Hz; high gamma, 100–200 Hz; and hfo; 200–400 Hz. For interregional analyses [phase-locking value (PLV), phase slope index (PSI), and Granger prediction], we focused on frequencies up to 100 Hz; spectral or time series data were subsequently downsampled to 250 Hz.
Anatomical reconstruction of recording sites.
Patients underwent preoperative, intraoperative, and postoperative imaging per routine clinical care. Preoperatively, stereotactic protocol MR images were obtained (Vario 3.0 T scanner, Siemens) that included T1- and T2-weighted sequences [T1, MPRAGE sequence: repetition time (TR), 2530 ms; echo time (TE), 2.85 ms; matrix size, 512 × 512 voxels; in-plane resolution, 0.5 × 0.5 mm2; 224 sagittal slices; 1 mm slice thickness; T2, SPACE sequence: TR, 3200 ms; TE, 409 ms; matrix size, 512 × 512 voxels; in-plane resolution, 0.5 × 0.5 mm2; 224 sagittal slices; 1 mm slice thickness]. Preoperative, intraoperative, and postoperative (in some cases) CT scans were also acquired [Lightspeed VCT Scanner for extraoperative imaging, GE Healthcare (tube voltage, 120 kV; tube current, 186 mA; data acquisition diameter, 320 mm; reconstruction diameter, 250 mm; matrix size, 512 × 512 voxels; in-plane resolution, 0.488 × 0.488 mm2; 267 axial slices; slice thickness, 0.625 mm); AIRO Mobile CT System Scanner for intraoperative scanning, Mobius Imaging (tube voltage, 120 kV; tube current, 240 mA; data acquisition diameter, 1331 mm; reconstruction diameter, 337 mm; matrix size, 512 × 512 voxels; in-plane resolution, 0.658 × 0.658 mm2; 182 axial slices; slice thickness, 1 mm)]. Postoperative MR images [Aera 1.5 T scanner, Siemens (T1, MPRAGE sequence: TR, 2300 ms; TE, 4.3 ms; matrix size, 256 × 256 voxels; in-plane resolution, 1.0 × 1.0 mm2; 183 axial slices; slice thickness, 1 mm; specific absorption rate, <0.1 W/g)] were typically obtained 1–2 d after the operation to confirm proper final electrode location.
To reconstruct recording locations, MR and CT images were coregistered using the WayPoint Navigator software (FHC). The raw DICOM images and the linear transform matrices were exported and applied to reconstructed image volumes using the AFNI command “3dAllineate,” bringing them into a common coordinate space (Cox, 1996; Li et al., 2016). Microelectrode depths were calculated by combining intraoperative recording depth information with electrode reconstructions obtained from postoperative images using methods described previously (Lauro et al., 2016, 2018). To determine the anatomic distribution of microelectrode recording sites across patients, preoperative T1-weighted MR images were registered to a T1-weighted MNI reference volume (MNI152_T1_2009c) using the AFNI command “3dQwarp” (Fonov et al., 2009). The resulting patient-specific transformation was then applied to recording site coordinates. MNI-warped recording coordinates were then assessed for proximity to structures such as the STN as delineated on the MNI PD25 atlas (Xiao et al., 2012, 2015, 2017). ECoG contacts were segmented from intraoperative CT volumes using the same DBStar processing as microelectrodes. Contacts were then projected onto individual cortical surface reconstructions generated from preoperative T1 volumes (Dale et al., 1999; Fischl et al., 2002; Saad and Reynolds, 2012; Trotta et al., 2018). Individual cortical surface reconstructions were coregistered to a standard Desikan–Destrieux surface parcellation (Argall et al., 2006; Desikan et al., 2006; Destrieux et al., 2010). Contacts were labeled and grouped as “premotor cortex” (PMC), “motor cortex” (MC), “somatosensory cortex” (SC), or “parietal cortex” if they contained the following anatomic parcellation labels: PMC, ctx_lh_G_front_sup, ctx_lh_G_front_middle; MC, ctx_lh_G_precentral; SC, ctx_lh_G_postcentral; and posterior parietal cortex (PPC), ctx_lh_G_parietal_sup, ctx_lh_G_pariet_inf-Supramar.
If a contact had more than one label (8 of 80 contacts), they were removed from further analysis.
Experimental design.
We used a visual–motor target tracking task to estimate the degree of motor dysfunction in a continuous fashion. Specifically, while patients with PD reclined on the operating bed in a “lawn-chair” position, a joystick was positioned within their dominant hand, and a boom-mounted display was positioned within their direct line of sight at a distance of ∼1 m. The task was implemented in MonkeyLogic (Asaad and Eskandar, 2008a,b; Asaad et al., 2013) and required subjects to follow a green target circle that moved smoothly around the screen by manipulating the joystick with the goal of keeping the white cursor within the circle (Fig. 1A). The target circle followed one of several possible paths (invisible to the subject), with each trial lasting 10–30 s. Each session consisted of up to 36 trials (∼13 min of tracking data), and subjects performed one to four sessions during the operation. Control subjects performed this task in an extraoperative setting.
Speed quantification.
To calculate movement speed, x- and y-joystick traces were 3 Hz low-pass filtered, and the euclidean change of cursor position was calculated over time. To standardize movement speed within patients, movement speed values within a session were minimum–maximum normalized into a measure of “slowness,” where 0 = highest speed and 1 = lowest speed.
Tremor amplitude quantification.
To calculate tremor, x- and y-joystick traces were 3–8 Hz bandpass filtered, and a one-dimensional linear projection of the filtered traces was calculated. Tremor amplitude and phase were calculated using the Hilbert transform of the resulting one-dimensional time series.
Tremor epoch design.
To standardize tremor amplitude across patients, tremor amplitude values from control subjects and patients were averaged into 4 s contiguous, nonoverlapping epochs. We chose our 4 s window size based in part on fMRI studies (Helmich et al., 2011), which based estimates of tremor amplitude on the timescale of echoplanar imaging TRs, which were 1–2 s. In addition, we calculated the autocorrelation of tremor amplitude within individual patient sessions and found that across all patients the central peak (>3 SDs above the mean) spanned 2 s. To capture the transition from one discrete state (no tremor) to another (sustained tremor), we chose a window size of 4 s to capture both states within one “tremor onset” epoch.
The resulting average and SD of the control tremor amplitude distribution were used to z-transform control subject and PD patient tremor amplitude epochs (Fig. 1B). To determine a cutoff to optimally differentiate control and PD population tremor data, receiver operating characteristic (ROC) tests were performed between supracutoff population data for cutoff values ranging from −2 (the lowest observed in both populations) to 10. The maximum area under the curve (AUC) value was observed for z = 3 (ROC AUC = 0.85), which was used for subsequent analyses.
“No tremor” and “sustained tremor” epochs were identified by 4 s epochs where the average tremor amplitude was subthreshold or suprathreshold. Potential tremor onset epochs were detected by taking the continuous tremor amplitude (1 ms samples) and identifying those time points where tremor amplitude crossed from subthreshold to suprathreshold levels. Epochs were then classified as “onset” if the mean of tremor amplitude samples in the subsequent 2 s were >3 SDs, and if the mean of tremor amplitude samples in the preceding 2 s was <3 SDs.
While there was no within-condition epoch overlap (i.e., each sustained tremor epoch was nonoverlapping), there was slight overlap between no tremor epochs with the pretrigger segment of tremor onset epochs (12 of 575 epochs across all subjects; mean ± SD of overlap, 1.275 ± 0.582 s). For sustained tremor, there were 18 of 171 epochs with some overlap with the post-trigger segment of tremor onset (mean ± SD of overlap, 1.008 ± 0.824 s).
Tremor frequency calculation and UPDRS correlation.
To calculate each patient's dominant tremor frequency (i.e., the frequency with the largest amplitude), a distribution of tremor amplitude was created by aggregating each patient's tremor amplitude epochs. In parallel, a frequency distribution was created by calculating the dominant (highest-power) tremor frequency within each epoch. A patient-specific dominant tremor frequency was then calculated as the frequency containing the highest aggregate tremor amplitude.
Correlations between task-derived tremor amplitude and UPDRS were conducted with subscores pertaining to the upper extremity relevant to the patient's task performance (rest tremor, postural tremor, finger taps, hand opening/closing, rapid alternating movements, rigidity). Each patient UPDRS subscore was Spearman correlated with the median of each patient's tremor amplitude distribution and was assessed for significance using a bootstrap null distribution (1000 iterations) where tremor medians were randomly shuffled with respect to UPDRS subscores.
Tremor/speed-spectral power correlation.
To determine whether spectral power across frequencies correlated with changes in tremor amplitude or slowness, linear mixed models were fit to 4 s epochs of averaged tremor/slowness and spectral magnitude of canonical frequency bands (theta, alpha, low beta, high beta, low gamma, high gamma, hfo). Models were fit within entire task sessions for each band and were specified as follows:
Tremor epoch spectral power modulation.
To determine whether the spectral band power at each structure differed by tremor epoch type, linear mixed models were used to compare spectral band power across epoch types by the following model:
Tremor-Neural theta phase-locking value.
To determine whether theta (3–8 Hz) in tremor and neural recordings were synchronized, the PLV was calculated with tremor and neural theta phase per trial (Lachaux et al., 1999). Theta-phase estimates for neural spectral data were calculated by taking the circular/angular mean for narrow-band phase estimates between 3 and 8 Hz at each time point (t), as follows:
To determine whether tremor-neural theta phase synchrony at each structure differed by tremor epoch type, linear mixed models were used to compare PLV values across epoch types by the following model:
As PLV and PPhC results were qualitatively similar, we reported PLV results.
Tremor–neural theta phase slope index.
To understand the lag–lead relationship between tremor (a bandpassed signal) and neural theta-phase locking, the PSI was calculated for the theta band (3–8 Hz) with 1 Hz frequency resolution (Nolte et al., 2008) using the “spectral_connectivity” Python toolbox (https://github.com/Eden-Kramer-Lab/spectral_connectivity, https://doi.org/10.5281/zenodo.4088934).
As the spectral_connectivity toolbox uses the multitaper transform for spectral analysis, the number of necessary tapers (L) was calculated by first calculating the time–half-bandwidth product (TW) using the desired frequency resolution (Δf; 1 Hz for parity with wavelet spectral analyses) and the time window of the entire trial (N, 4 s; Prerau et al., 2017), as follows:
We subsequently used TW to calculate L using the floor function (
With our parameters, three Slepian tapers were used for whole-trial single-window PSI estimates, as follows:
PSI was then estimated from the imaginary
Phase offsets between 1 Hz frequency bands (Δf) within theta (F) were used to calculate the phase slope. Because of our short-timescale windowed application of PSI, we did not normalize values of PSI by their SD (Young et al., 2017). To determine whether tremor or neural recordings exhibited directional theta influence, the empirical PSI was compared with a null distribution of 1000 PSI values generated from shuffling time series of one signal across trials. The p values were calculated empirically from the resulting distribution and corrected for multiple comparisons with the Benjamini–Hochberg method at q = 0.05.
Tremor epoch interregional phase-locking value.
To compare time-varying phase synchrony across structures, the PLV was calculated across each structure pair (j, k) per 1 Hz frequency band (f) from 1 to 100 Hz using wavelet-derived spectral data, as follows:
To determine whether pairwise frequency band PLVs differed by tremor epoch type, linear mixed models were used to compare PLV values across epoch types by the following model:
Tremor epoch interregional Granger prediction.
To understand whether tremor epoch-related dynamic changes in spectral power or synchrony were driven by dynamic directional influences of one structure onto another, nonparametric spectral Granger prediction (GP) was calculated between each structure pair using the “spectral_connectivity” python toolbox. Specifically, frequency information (1 Hz frequency resolution) for each structure–time series pair was calculated using a single 4000 ms multitaper window (three tapers). From there, a frequency-based estimation of information flow between structures was calculated using a cross-density spectral matrix (Dhamala et al., 2008). Subsequently, frequency-specific (f) GP (i.e., the log-ratio of total frequency power over nonpredicted frequency power) was calculated between structure pairs (j, k) for each epoch type using S, the spectral transfer matrix (H), and the noise covariance matrix (∑), as follows:
To determine whether one structure exhibited frequency-specific Granger prediction on another, the empirical GP was compared with a null distribution of 1000 GP values generated from shuffling time series of one structure across trials. The p values for each frequency were calculated empirically from the resulting distribution and corrected for multiple comparisons with the Benjamini–Hochberg method at q = 0.05.
To understand how GP varied as a function of time, frequency information for each structure–time series pair were calculated in 2000 ms windows with 100 ms overlap using the multitaper transform for each event trial. To maintain the same number of tapers (three tapers) between static and dynamic GP analyses, frequency resolution was increased to 2 Hz for dynamic GP calculation. To determine whether one structure exhibited time-varying directional influence on another, the empirical GP was compared with a null distribution of 1000 GP values generated from shuffling time series of one structure across trials. The p values for each time and frequency point were calculated empirically from the resulting distribution and corrected for multiple comparisons with the Benjamini–Hochberg method at q = 0.05. Resulting significant time–frequency clusters were additionally filtered by considering only clusters whose area was greater than the 95th percentile of all Benjamini–Hochberg method-corrected significant clusters.
Tremor epoch interregional phase slope index.
To calculate theta directed connectivity across structures, the PSI was used for the theta band (3–8 Hz) with 1 Hz frequency resolution across structures. Frequency information (1 Hz frequency resolution) for each structure–time series pair was calculated in a single 4000 ms window using the multitaper transform (three tapers) for each event trial. To determine whether one structure exhibited PSI influence on another, the empirical PSI was compared with a null distribution of 1000 PSI values generated from shuffling one structure's time series across trials. The p values were calculated empirically from the resulting distribution and corrected for multiple comparisons with the Benjamini–Hochberg method at q = 0.05.
To calculate the time-varying PSI between broad frequency bands, PSI was calculated using a 2000 ms window sliding by 100 ms (three tapers with 2 Hz frequency resolution). A bootstrap was then performed, and empirical p values for each time point were corrected for multiple comparisons with the Benjamini–Hochberg method at q = 0.05.
Statistical analysis.
Data in text are represented as mean ± SD. Because data were aggregated across multiple subjects, linear mixed models performed with the “statsmodels” Python toolbox were used to disentangle the fixed effects of subject population, event condition, or spectral band power from the random effects of each subject's dataset (Lindstrom and Bates, 1988; Seabold and Perktold, 2010). All linear mixed models were random intercept models, where each subject's dataset was assigned a random intercept 1|Subject. Once a model was fit, p values were calculated from z-scored parameter estimates (parameters estimates divided by their SEs) against the normal distribution. Because directed connectivity measures (PSI, GP) use multiple epochs for a single estimate of directed connectivity, linear mixed models were not able to account for individual patient variability in these results. Instead, we used bootstrapping where recordings were shuffled across all epochs aggregated across all subjects, and p values were calculated empirically from the resulting distribution. All other statistical tests, unless otherwise specified, were conducted in the “scipy” Python environment. The p values were controlled for multiple comparisons by using the Benjamini–Hochberg procedure at q = 0.05 (Benjamini and Hochberg, 1995).
Data availability.
The datasets supporting the current study have not been deposited in a public repository because they contain patient information but are available along with analysis code on request.
Results
Intraoperative behavioral and neural data acquisition
Ten patients with PD undergoing DBS implantation and 14 age-matched control subjects (see Materials and Methods) performed a simple visual–motor task where they followed an onscreen target using a joystick-controlled cursor with their dominant hand (Fig. 1A). Each patient performed one to four sessions of this target-tracking task during the procedure for a total of 27 sessions, while control subjects performed 1 session each for a total of 14 sessions. Tremor amplitude and cursor speed were quantified continuously from the x- and y-joystick traces (control subjects, n = 1856 epochs; PD, n = 3400 epochs). As patients were not all clinically tremor dominant, and because all subjects contributed variable amounts of data, linear mixed models were used to quantify the difference of tremor/speed across subject populations while accounting for individual subject variability. While the resulting PD and control speed distributions were distinct (linear mixed-model coefficient = 0.884, z = – 3.138, p = 0.002), PD distributions trended toward having increased tremor relative to control subjects (linear mixed-model coefficient = 0.180, z = 1.859, p = 0.063; Fig. 1B,C). Nevertheless, the partial overlap of the PD and control tremor distributions (indicative of periods without tremor in PD patients), along with the long right tail of the PD distribution, gave us a large dynamic range of tremor to analyze with respect to neural data. The dominant tremor frequency across patients was 4.48 ± 0.57 Hz. While tremor amplitude correlated with the resting tremor UPDRS subscore across patients (Spearman ρ = 0.92, p < 0.001, bootstrap test), it did not with the postural tremor subscore (ρ = 0.54, p = 0.065, bootstrap test). Based on the distinct tremor frequency peak and its correlation with clinical measures of resting tremor, we interpreted our task-derived tremor as reflective of resting tremor (Dirkx et al., 2018).
Tremor and movement speed calculated from the intraoperative visual–motor task. A, Left, Schematic of task target (green) and joystick (gray) traces from a single trial. Middle top, Bandpass-filtered x- and y-joystick traces from the task trial. Middle bottom, Low-pass-filtered x- and y-joystick traces from the task trial. Right top, One-dimensional projection of bandpass-filtered traces (black), with tremor amplitude measured from the envelope (orange). Right bottom, Cursor speed measured from low-pass-filtered traces (black). B, Distribution of 4 s tremor amplitude epochs for control subject and PD patient populations. ∘, Degrees of visual angle. Vertical dashed line indicates ROC-derived cutoff value between control and PD populations. While there is overlap on the left side of the distribution (patients with PD can exhibit control-like performance), the PD distribution is highly skewed on the right side of the distribution, allowing a large range of tremor expression. C, Distribution of 4 s speed epochs for control subject and PD patient populations. The bimodality of the control distribution corresponded to the preprogrammed speed of the onscreen target. Despite this, note that the PD distribution is shifted toward lower speed values. D, Coronal view of microelectrode recording density on an MNI reference volume. The inset panel displays a close-up view of the subthalamic nucleus (outlined in black). L, Left. E, Recording density of ECoG contacts on an MNI reference surface.
Across the 10 patients with PD, we obtained 81 microelectrode recordings within the STN (peak recording density: MNI coordinates,
Tremor is a neurophysiologically distinct motor feature of Parkinson's disease
To understand the relationship of broadband neural activity to tremor expression, we examined the correlation between tremor amplitude and spectral power in neurophysiological recordings. Sorting session-wide spectral data by tremor epochs (rather than according to time) revealed informative band-specific patterns of activity (STN, n = 81 session recordings; PMC, n = 75; MC, n = 49; SC, n = 51; PPC, n = 41; Fig. 2A). Specifically, across cortical structures with the exception of PMC, spectral power in low and high beta frequency ranges (12–30 Hz) were found to negatively correlate with tremor amplitude (linear mixed-model coefficients = −0.325 to −0.902, z = −5.000 to −18.931, p ≤ 5.77 * 10–7; Fig. 2B). Interestingly, beta power appeared to drop off fairly quickly with even low levels of tremor becoming evident (SC: power curve fit, r2 = 0.77; linear fit, r2 = 0.54). Meanwhile, theta power positively correlated with tremor amplitude in the STN, MC, and SC (linear mixed-model coefficients = 0.076–0.732, z = 3.875–6.569, p ≤ 1.07 * 10–4).
Tremor and slowness exhibit distinct spectral power correlations with intracranial recordings. A, Population-averaged task session spectral power, sorted by the tremor amplitude of each epoch (left) or slowness (right). For ease of visualization, frequency power was z-scored within frequencies across epochs. B, Average session-wide narrowband (1 Hz) spectral Spearman correlation (ρ) with tremor amplitude and slowness. Note that while beta frequencies exhibited an opposing relationship with tremor and slowness, theta frequencies exhibited a distinct positive correlation with tremor.
To compare tremor-related neural activity with a distinct PD motor feature (specifically bradykinesia), neural data were also analyzed with respect to movement “slowness” during the same target-tracking task. Note that PD subjects appeared to lack a higher mode of movement velocity that was clearly present in control subjects, reflecting an inability to move the cursor consistently as quickly as the target (Fig. 1C). We calculated a minimum-maximum normalized measure of inverse cursor speed (0 = highest speed, 1 = lowest speed) to capture this effect as a positive pathologic sign, parallel to the sign of tremor. In contrast to tremor, we observed positive correlations between slowness and alpha/beta (8–30 Hz) power in all cortical structures (linear mixed-model coefficients = 0.159–1.141, z = 6.937–20.587, p ≤ 4.01 * 10–12; Fig. 2B). However, theta did not show a significant correlation with slowness in any structure (p > 0.05). Thus, theta appeared to relate specifically to tremor, whereas the relationship to beta activity was generally reversed between these PD-related motor manifestations. So while there was broadly the appearance of a symmetric opposition between tremor and slowness in terms of their correlations with neural activity across frequencies (Fig. 2B), this difference in the theta frequency relationship, as well as perhaps a consistent difference in mid-gamma (in which the correlation with tremor was typically close to 0, but the correlation with slowness was typically greater in magnitude and negative in direction), suggest that these motor features are not simply opposite ends of a single spectrum but rather have distinct fingerprints in neural activity.
Subthalamic theta preceded tremor at onset
Because lower-frequency oscillations, particularly theta, were most consistently and strongly positively associated with tremor across structures, and because they encompassed the range of observed tremor frequencies from a behavioral perspective (4–6 Hz), we next turned our attention to understanding the relationship of theta band activity within each structure to tremor-defined epochs. Using a control versus PD subject ROC-derived tremor threshold (see Materials and Methods), behavioral and spectral data were organized into 4 s epochs and categorized as follows: no tremor epochs (n = 575 epochs, 2300 s), tremor-onset epochs (n = 406 epochs, 1624 s), and sustained tremor epochs (n = 171 epochs, 684 s; Fig. 3A). All 10 patients contributed at least one epoch to the No Tremor and Tremor Onset conditions, with 6 patients contributing at least one epoch to the Sustained Tremor condition. The resulting behavioral and spectral data were aggregated across subjects (No Tremor: n = 1725 tremor–STN paired recording epochs, Tremor Onset: n = 1218, Sustained Tremor: n = 513).
Spectral power during different tremor dynamic states. A, Tremor event design. Based on a population-based tremor ROC threshold, epochs representing different states of tremor dynamics were isolated. For each event type, the average tremor amplitude (±SE) in patients with PD relative to control subjects is displayed over time. Horizontal dashed line denotes the tremor threshold (3 SDs relative to control subjects). Vertical dashed line (t = 0) in tremor-onset events represents the “trigger” where tremor amplitude crossed the tremor threshold. B, Average spectral power (±SE) across frequencies for each tremor event type, by recording site. Vertical dashed lines represent frequency band borders. While theta oscillations increased in power across STN, MC, and SC, increased tremor was associated with increased alpha/low beta power in PMC and PPC.
STN theta power was indeed significantly elevated during tremor onset (linear mixed-model coefficient = 0.070, z = 8.039,
Neural theta exhibited structure-specific temporal relationships with tremor. A, Histograms of per-trial PLVs between tremor and neural theta by tremor state. Solid lines indicate normal distribution fit to each tremor state PLV histogram, while vertical dashed lines indicate the median of each tremor state PLV histogram. The y-axis indicates the proportion of trials within each PLV histogram bin. Note that STN histograms for tremor onset and sustained tremor are highly overlapping. B, PSI between tremor and neural theta by tremor state. Positive values indicated that tremor phase preceded neural phase, while negative values indicated that neural phase preceded tremor. Magenta asterisks indicate significant (p < 0.05, bootstrap test) PSI effects.
In light of this close relationship between STN theta and tremor, we next examined the temporal relationship between STN theta and tremor phase. Specifically, we calculated the PSI between tremor and STN theta phase. Because the PSI considers multiple phase relationships within a range of frequencies, it can succeed in determining the net leading or lagging oscillation in a manner that avoids the circularity problem inherent in methods such as the PLV (Nolte et al., 2008). Here, the PSI revealed STN theta led tremor exclusively during tremor onset (p = 0.011, bootstrap test; Fig. 4B), consistent with a causal role for the STN in the initiation but not necessarily the maintenance of tremor.
Somatosensory cortex theta consistently followed tremor
Like the STN, SC theta power positively correlated with tremor amplitude. Therefore, we investigated whether this spectral–tremor relationship varied similarly with tremor state (No Tremor, n = 1256 tremor–SC paired recording epochs; Tremor Onset, n = 746; Sustained Tremor, n = 150). SC theta power was indeed significantly elevated during tremor onset (linear mixed-model coefficient = 0.010, z = 4.831,
However, in contrast to the STN, phase-slope analysis of tremor and SC theta phase revealed that SC theta phase followed tremor phase during both tremor onset and sustained tremor (p ≤ 0.002, bootstrap test; Fig. 4B). Therefore, the strong tremor-related theta oscillation seen in SC was reflective rather than causal of tremor.
Motor cortex theta consistently preceded tremor
Like the STN and SC, MC theta power showed a clear graded relationship with tremor magnitude (Fig. 2). Examining MC theta power across tremor states (No Tremor, n = 1066 tremor–MC paired recording epochs; Tremor Onset, n = 692; Sustained Tremor, n = 312) revealed that it was relatively increased during tremor onset (linear mixed-model coefficient = 0.006, z = 2.701, p = 0.007) but not sustained tremor (linear mixed-model coefficient = 0.002, z = 0.429, p = 0.668; Fig. 3B). Furthermore, MC-tremor theta PLV increased from no-tremor to tremor-onset (linear mixed-model coefficient = 0.018, z = 2.646, p = 0.008) to sustained tremor (linear mixed-model coefficient = 0.105, z = 10.184,
Tremor-related theta transitioned from STN to cortex during tremor onset
Because both STN and MC theta power were elevated during tremor onset, and STN and MC theta phase led tremor phase during tremor onset, we investigated the dynamics of STN–MC coupling during the dynamics of tremor initiation (No Tremor, n = 3198 STN–MC paired recording epochs; Tremor Onset, n = 2076; Sustained Tremor, n = 936). Static phase slope analysis of STN and MC revealed that STN theta led MC theta during tremor onset (p < 0.001, bootstrap test; Fig. 5A). To understand whether this phase relationship was time locked to increasing tremor, we calculated STN-MC theta PSI as a function of time within the tremor onset window. Within this epoch, STN theta preceded MC theta most consistently ∼0.5 s after tremor detection (t = 0) to the end of the tremor onset epoch (t = 0.5–1.0 s; p < 0.05, bootstrap test; Fig. 5B). At no point in this window did MC theta appear to precede STN theta.
Tremor initiation was driven by the subthalamic nucleus. A, Static PSI between STN and MC recordings during tremor states. Magenta asterisks indicate significant (p < 0.05, bootstrap test) PSI effects. B, Dynamic PSI between STN and MC theta during tremor onset. Highlighted regions indicate significant PSI (p < 0.05, bootstrap test). Vertical dashed line (t = 0) indicates tremor-onset trigger. C, Directed GP between STN and MC theta during tremor onset. Vertical dashed line (t = 0) indicates tremor-onset trigger. Highlighted regions indicate significant Granger prediction (p < 0.001, bootstrap test). D, Static PSI between STN and SC recordings during tremor states. Magenta asterisks indicate significant (p < 0.05, bootstrap test) PSI effects. E, Dynamic PSI between STN and SC theta during tremor onset. Highlighted regions indicate significant PSI (p < 0.05, bootstrap test). Vertical dashed line (t = 0) indicates tremor-onset trigger. F, Directed GP between STN and SC theta during tremor onset. Vertical dashed line (t = 0) indicates tremor-onset trigger. Highlighted regions indicate significant Granger prediction (p < 0.001, bootstrap test). G, Summary of theta PSI results. Solid lines represent directed functional connectivity between neural regions and tremor.
We also investigated whether STN theta and MC theta power influenced each other by calculating time-varying nonparametric spectral GP (see Materials and Methods). Briefly, a nonzero GP at a particular frequency indicated that spectral power in one structure was predictive of spectral power in another. Unlike the PSI, GP allows the disentangling of asymmetric, bidirectional influences across two signals (Dhamala et al., 2008). As with PSI, STN theta power predicted MC theta power from 200 ms after the tremor onset trigger to the end of epoch (t = 0.2–1.0 s; p < 0.05, bootstrap test; Fig. 5C). Again, MC theta did not predict STN theta at any point in the epoch. Together, these results converged to suggest STN theta drove MC theta during tremor onset.
Once tremor was established, however, the theta phase slope relationship flipped, with MC theta phase preceding STN theta phase (Fig. 5A), revealing a dynamic transition with increasing tremor. Together with the loss of STN theta influence over tremor during sustained tremor (Fig. 4B), tremor output appeared to become cortically driven rather than STN driven as tremor became established.
Because the STN and SC both exhibited positive correlations between theta power and increasing tremor, we also investigated whether STN/SC dynamics varied during tremor onset (No Tremor, n = 3768 STN–SC paired recording epochs; Tremor Onset, n = 2238; Sustained Tremor, n = 450). Like MC, static phase slope analysis of STN and SC theta revealed that STN theta led SC during tremor onset (p < 0.001, bootstrap test; Fig. 5D). Dynamic STN–SC PSI additionally revealed that STN theta led SC theta between 200 ms after the tremor onset trigger to the end of the epoch (t = 0.2–1.0 s; p < 0.001, bootstrap test; Fig. 5E). Simultaneously, STN theta power predicted SC theta power from 400 ms before the tremor-onset trigger to end of the tremor-onset epoch (t = – 0.4 to +1.0 s; p < 0.001, bootstrap test; Fig. 5F). During sustained tremor epochs, however, the theta-phase slope relationship between STN and SC became ambiguous (p = 0.091, bootstrap test), again representing a loss of STN influence over cortical theta activity (Fig. 5D). Altogether, although the STN drove both tremor and cortical theta as tremor emerged, the transition to sustained tremor was accompanied by a decoupling of the STN from cortex in the theta band (Fig. 5G).
Motor cortex lost influence over posterior cortices with increasing tremor
As STN–MC theta phase influence flipped from tremor onset to sustained tremor, we investigated whether the functional connectivity of MC extended to other cortical regions with increasing tremor. To understand whether tremor-mediated corticocortical interactions occurred in frequency bands other than theta, we calculated both nondirected (PLV) and directed (GP) functional connectivity between the MC and other cortical regions across the 3–100 Hz spectrum (paired recording epochs from No Tremor, Tremor Onset, and Sustained Tremor conditions: MC–PMC: n = 2692, 2064, 757; MC–SC: n = 2190, 1130, 210; MC–PPC: n = 1458, 1074, 935). MC–SC PLV across any frequency band did not modulate by tremor state (PLV; linear mixed model, p > 0.05; Fig. 6A). To identify whether synchrony detected by the PLV was driven by one structure in the pair, broad-spectrum GP was calculated. In the absence of tremor, we found that MC predicted SC high beta/low gamma power (p < 0.001, bootstrap test) and SC predicted MC theta/alpha/low beta/mid-gamma power (p < 0.001, bootstrap test; Fig. 6B). During sustained tremor, however, MC theta now predicted SC theta (GP, 2.34-fold difference, p < 0.001, bootstrap test), and SC low beta/low gamma predicted MC low beta/low gamma (GP, 1.99- to 2.18-fold difference; p < 0.001, bootstrap test).
During sustained tremor, gamma coupling between premotor/motor and somatosensory/parietal cortices decreased. A, PLV between MC and other cortical regions. Lines ± shaded borders represent average ± SE PLV. Highlighted frequency ranges indicate increased (orange) or decreased (blue) PLV with increasing tremor. B, Pairwise GP between MC and other cortical regions. The title of each subpanel indicates the directionality of the structure pair GP. Highlighted frequency ranges indicate increased (orange) or decreased (blue) GP with increasing tremor. Note that MC broad-spectrum coupling with SC and PPC generally decreased with increasing tremor. C, PLV between SC and other cortical regions. Lines ± shaded borders represent average ± SE PLV. Highlighted frequency ranges indicate increased (orange) or decreased (blue) PLV with increasing tremor. D, Pairwise GP between SC and other cortical regions. The title of each subpanel indicates the directionality of the structure pair GP. Highlighted frequency ranges indicate increased (orange) or decreased (blue) GP with increasing tremor. Note that tremor generally shifted the frequency of coupling among SC, PPC, and PMC from gamma to alpha/low beta with increasing tremor. E, PLV between the STN and cortical regions. Lines ± shaded borders represent average ± SE PLV. Highlighted frequency ranges indicate increased (orange) PLV with increasing tremor. F, Pairwise GP between the STN and cortical regions. The title of each subpanel indicates the directionality of the structure pair GP. Highlighted frequency ranges indicate increased (orange) GP with increasing tremor, and increased GP specific to Tremor Onset (gray). Note that while broad-spectrum STN–cortical PLV generally increased with increasing tremor, directional changes were less distinct. For ease of visualization, GP curves were low-pass filtered and frequencies within 58–62 Hz were masked. Vertical dashed lines represent frequency band borders.
MC–PPC PLV similarly did not modulate as tremor increased (PLV linear mixed model, p > 0.05; Fig. 6A). However, Granger analysis revealed that PPC theta/alpha/low beta predicted MC theta/alpha/low beta regardless of tremor state (GP, 1.02–4.20-fold difference; p < 0.001, bootstrap test in all tremor states; Fig. 6B). In contrast, while MC high beta/low gamma/mid-gamma predicted PPC high beta/low gamma/mid-gamma in the absence of tremor (GP, 1.07–1.48-fold difference; p < 0.001, bootstrap test), this relationship flipped during sustained tremor, with PPC high beta/low gamma/mid-gamma predicting MC high beta/low gamma/mid-gamma (GP, 1.50–1.99-fold difference; p < 0.001, bootstrap test).
In sum, MC exerted less influence over posterior (SC, PPC) and anterior (PMC) cortical regions with increasing tremor. Specifically, MC–PMC PLV decreased within low beta/low gamma (12–20/30–60 Hz) specifically during sustained tremor (PLV linear mixed-model coefficients, −0.010 to −0.017; z = −2.100 to −2.323, p ≤ 0.035; Fig. 6A). While not within the same frequency range, PMC theta/alpha also appeared to predict MC theta/alpha during sustained tremor (GP, 2.68- to 2.71-fold increase; p < 0.001, bootstrap test; Fig. 6B).
Premotor cortex coupled with posterior cortices during tremor
Because SC decoupled from the STN during sustained tremor while still reflecting tremor output, we investigated whether SC instead coupled with other cortical regions as tremor increased (paired recording epochs from No Tremor, Tremor Onset, and Sustained Tremor conditions: SC–PMC: n = 3442, 2292, 459; SC–PPC: n = 1284, 808, 165; PMC–PPC: n = 1780, 1462, 1029). SC–PPC PLV similarly did not modulate as tremor increased (PLV, linear mixed model, p > 0.05; Fig. 6C). Although PLV did not significantly modulate with tremor, PPC theta/alpha/low beta (8–20 Hz) predicted SC theta/alpha/low beta during sustained tremor (GP, 1.61- to 8.55-fold difference; p < 0.001, bootstrap test; Fig. 6D). While SC low gamma/mid-gamma predicted PPC low gamma/mid-gamma during the absence of tremor (GP, 1.23- to 1.32-fold difference; p < 0.001, bootstrap test), this did not hold for sustained tremor. Thus, SC–PPC connectivity shifted to a distinct state during sustained tremor, with PPC predicting lower frequencies (theta, alpha, low beta) in SC. At the same time, higher frequency (gamma) directed connectivity between SC and PPC decreased as tremor increased.
SC and PMC interactions exhibited decreases in functional connectivity, with decreased high beta–low gamma (20–60 Hz) PLV (PLV linear mixed-model coefficients, −0.022 to −0.007; z = −1.975 to −4.488; p < = 0.048) with increasing tremor (Fig. 6C). Like PPC, SC alpha/low beta was driven by PMC alpha/low beta specifically during sustained tremor (GP, 4.41- to 13.10-fold difference; p < 0.001, bootstrap test; Fig. 6D). Thus, in contrast to MC, which lost influence over posterior cortical regions, SC became increasingly influenced by both posterior (PPC) and anterior (PMC) cortices with increasing tremor. However, this increase in connectivity was specific to alpha/low beta frequencies while gamma coupling decreased between SC and PMC/PPC.
To follow the spread of tremor-related cortical coupling, we investigated whether PMC and PPC interacted during sustained tremor. Here, we observed an exaggerated version of the same tremor-induced frequency shift (gamma to beta) of power and phase synchrony. When analyzing tremor epoch-related spectral power in PMC and PPC in Figure 3B, both regions demonstrated tremor-related decreases in low beta/high beta frequencies (linear mixed-model coefficients, −0.007 to −0.019; z = −2.358 to −5.256; p ≤ 0.018). At the same time, PMC exhibited decreases in low gamma power during sustained tremor relative to no tremor (linear mixed-model coefficients, −0.003; z = – 2.895; p = 0.004).
These similar changes in power were mirrored by changes in PMC–PPC PLV synchrony (Fig. 6C). PMC–PPC low to mid-gamma PLV decreased as tremor increased (PLV linear mixed-model coefficients, −0.016 to−0.020; z = −3.021 to −3.367; p ≤ 0.003), while PMC–PPC alpha/low beta PLV increased with tremor (PLV linear mixed-model coefficients, 0.018–0.020; z = 2.348–3.253; p ≤ 0.018). Regardless of tremor state, PMC–PPC phase synchrony was driven by PMC onto PPC. When tremor was absent, PMC low gamma/mid-gamma predicted PPC low gamma/mid gamma (GP, 1.98- to 2.12-fold difference; p < 0.001, bootstrap test; Fig. 6D). During sustained tremor, PMC low beta/high beta power predicted PPC low beta/high beta power (GP, 2.87- to 6.26-fold difference; p < 0.001, bootstrap test).
Overall, tremor was associated with a frequency shift (gamma to beta) of power and phase synchrony among PMC, PPC, and SC. Specifically, PMC exerted increasing influence over posterior regions (SC, PPC) in lower frequencies (alpha, low beta) with increasing tremor. However, this increase in lower-frequency coupling coincided with decreases in higher-frequency coupling (gamma). In addition, directional gamma influence between MC and PPC flipped with increasing tremor (MC → PPC in the absence of tremor, PPC → MC during sustained tremor), revealing that sustained tremor is a state of altered gamma synchrony across sensorimotor cortex.
STN broadly synchronized with, but selectively influenced, sensorimotor cortex during sustained tremor
Finally, to understand whether STN influence over cortex extended beyond theta, functional and directed connectivity were calculated between the STN and sensorimotor cortex (paired recording epochs from No Tremor, Tremor Onset, and Sustained Tremor conditions: STN–PMC: n = 4680, 3732, 1281; STN–MC: n = 3198, 2076, 936; STN–SC: n = 3768, 2238, 450; STN–PPC: n = 2154, 1698, 1437). STN–cortical theta PLV synchrony (with the exception of SC) increased as a function of tremor (PLV linear mixed-model coefficients, 0.012–0.025; z = 2.873–6.827; p ≤ 0.004; Fig. 6E). However, directed connectivity between the STN and cortex was specific to tremor state (Fig. 6F). STN theta power (4–6 Hz) predicted both MC theta power (GP, 2.29-fold difference; p < 0.001, bootstrap test) and SC theta power (GP, 1.79-fold difference; p < 0.001, bootstrap test) exclusively during tremor onset. In contrast, STN theta power predicted PPC theta power (GP, 1.86-fold difference; p < 0.001, bootstrap test) and PMC theta power (GP, 1.59-fold difference; p < 0.001, bootstrap test) only during sustained tremor. Thus, consistent with the PSI results, the STN shifted its influence over cortex in the theta band (STN → MC/SC during tremor onset; STN → PMC/PPC during sustained tremor) across dynamic tremor states.
Discussion
Using a naturalistic behavioral task, we were able to characterize tremor dynamics and to isolate specific tremor states, particularly tremor onset and maintenance. Across structures, we found that theta power positively and beta power negatively correlated with tremor, as has been found in previous reports (Hirschmann et al., 2013; Qasim et al., 2016; Asch et al., 2020). However, our study is the first to dissect electrophysiological correlates of tremor onset and sustained tremor. During the emergence of tremor, not only did STN and motor cortical theta power increase, but STN and motor cortical theta phase preceded the phase of tremor. Moreover, STN theta activity drove motor cortical theta during tremor onset, suggesting a direct role of the STN in initiating tremor output.
Once tremor emerged, however, motor cortex appeared to sustain tremor. At the same time, motor cortex became less coupled with somatosensory and parietal cortices, despite the presence of prominent somatosensory cortex theta power, which closely followed tremor. Instead, premotor cortex synchronized via low beta frequencies with posterior cortices (somatosensory, parietal) at the expense of gamma frequency synchronization observed in the absence of tremor. This low beta synchrony was notably asymmetric across these structures, with premotor cortex exerting influence over posterior cortices.
Together, although tremor amplitude corresponded to global changes in theta and beta power, the relationship between these frequency bands to tremor output was highly structure specific. While STN–motor cortical interactions appeared to initiate tremor, premotor cortex-driven network effects may help sustain tremor. This STN-mediated dynamic reorganization of cortical connectivity is consistent with both the “dimmer switch” model and the intrinsic and extrinsic cortical loops of Parkinson's tremor (Helmich et al., 2011; Volkmann et al., 1996; Fig. 7). Like the GPi, we revealed that the STN acted as a “switch” to mediate the onset of tremor by influencing motor cortex (Dirkx et al., 2016). While these STN–motor cortical interactions formed the intrinsic loop of tremor output, we expanded this model to reveal that shifts from gamma to beta synchrony across premotor–parietal cortices reflected the extrinsic loop in the stable tremor state.
Synthetic model of subcortical–cortical interactions during tremor. Solid lines represent directed functional connectivity between neural regions and tremor. Dashed lines during Sustained Tremor represent interactions from the No Tremor state that are no longer present.
Tremor onset was mediated by subthalamic theta driving motor cortex
STN theta amplitude positively correlated with tremor amplitude regardless of tremor dynamic states. While the phase of STN theta consistently preceded tremor phase during tremor onset, it did not during sustained tremor. However, STN theta activity was still significantly phase locked to tremor during sustained tremor. This mixed relationship to tremor may reflect several roles of STN: interconnections with GPi contribute to tremor initiation, while disynaptic connections with cerebellum may influence ongoing monitoring of tremor output (Bostan et al., 2010; Helmich et al., 2011). Indeed, STN projections to cerebellar cortex may perhaps propagate tremor–frequency oscillations within the basal ganglia to motor cortical–thalamo–cerebellar loops (Wu and Hallett, 2013; Bostan and Strick, 2018). Future experiments combining STN and cerebellar recordings could describe this tremor-onset mechanism, while trying to disentangle the neural control of tremor amplitude and phase (Cagnan et al., 2014; Helmich et al., 2021).
Regardless, STN theta drove motor cortex activity during tremor onset. Although tremor has previously been found to decrease beta coherence between STN and motor cortex (Qasim et al., 2016), while increasing theta coherence (Hirschmann et al., 2013), we demonstrated directed theta-phase interactions from STN to motor cortex specifically during tremor onset. While a previous case study of tremor onset displayed local STN and cortical alpha/beta power changes with tremor onset (Hirschmann et al., 2019), we show here that STN and motor cortical theta activity are directionally linked. We also demonstrated that during sustained tremor, the STN–motor cortex theta phase slope relationship reversed, suggesting the theta influence over sustained tremor shifted source from STN to cortex.
Motor cortex desynchronized with posterior cortices while sustaining tremor
As tremor progressed, motor cortex theta increasingly drove tremor. While previous studies have correlated motor cortical activity to tremor (Timmermann et al., 2003; Helmich et al.,2011), this is the first study, to our knowledge, that has demonstrated a directed relationship between ECoG recordings and tremor. Although motor cortex was synchronized to tremor, motor cortex appeared to desynchronize with other cortical structures with the exception of premotor cortex, as has been found previously (Timmermann et al., 2003; Qasim et al., 2016). While other studies have found that motor cortex increased its synchrony with premotor and parietal cortices during tremor (Hirschmann et al., 2013), this was calculated only at tremor and double-tremor frequencies.
Tremor reorganized premotor and parietal cortical coupling
Although premotor and parietal cortices did not exhibit a direct theta relationship to tremor, changes in tremor initiated a frequency shift in premotor–parietal coupling dynamics. In the absence of tremor, these regions were functionally coupled at higher frequencies (high beta, low to mid-gamma). fMRI studies in patients with PD have found that these regions exhibit overactive BOLD activity during self-initiated sequential hand movements (Samuel et al., 1997), which is hypothesized to compensate for decreased BOLD activity in frontostriatal circuits in the dopamine-depleted state (Wu et al., 2011). Furthermore, cortical gamma frequency power and synchrony are associated specifically with voluntary movement (Crone et al., 1998; Miller et al., 2007). In our study, this bidirectional premotor–parietal gamma activity may have reflected task monitoring and spatial tracking (motor output) using sensory information.
During sustained tremor, however, parietal and premotor cortices both exhibited increases in low beta power. This low beta activity was also functionally coupled, with premotor driving parietal cortex. Elevated low beta oscillations have been observed in premotor cortex recordings in MPTP nonhuman primates with predominantly akinetic/rigid symptoms (Wang et al., 2017). While not observed in our study, increased premotor high beta influence over the STN has also been found to correlate with akinetic/rigid symptoms (Sharott et al., 2018). Premotor low beta oscillations may function here in a similar antikinetic fashion with other cortical structures during tremor.
In any case, with increasing tremor premotor–parietal gamma activity diminished while premotor low beta activity drove parietal activity. These frequency shifts may be best understood in the framework of communication-through-coherence theory (Fries, 2015). Specifically, while symmetric or bottom-up gamma oscillations permit effective and precise transmission of motor-related information across structures, lower-frequency oscillations such as alpha/beta act as top-down feedback. Here, task-related gamma synchrony observed across sensorimotor cortex decreased with tremor. In contrast, lower-frequency oscillations such as low beta increased in synchrony, perhaps reflecting an absence of voluntary movement, which normally acts to suppress tremor (Naros et al., 2018).
Implications for closed-loop deep brain stimulation
Because of the clinical interest in developing adaptive closed-loop DBS to more precisely treat PD symptoms, various electrophysiological observations have been investigated as potential tremor biomarkers to inform stimulation (Hirschmann et al., 2017; Shah et al., 2018; Yao et al., 2020). While promising, the features used for tremor detection do not take into account the dynamic nature of tremor—namely, the distinct neurophysiological signature of tremor onset. Because of the breadth of STN beta-frequency oscillation research in PD, initial closed-loop DBS efforts have focused on using beta oscillations as a proxy for bradykinesia symptoms (Little et al., 2013, 2016a,b; Velisar et al., 2019). However, beta-driven DBS has been shown to worsen tremor in some patients (He et al., 2020; Piña-Fuentes et al., 2020).
Here, we demonstrated that subthalamic theta was present whether tremor was emerging or sustained. The addition of STN theta-based biomarkers to closed-loop DBS could help treat the separate symptom axis of tremor. Further, we have provided the best evidence to date that cortical ECoG theta is a robust marker for tremor. Specifically, we found that motor cortical theta was synchronized to STN theta during tremor states, and that somatosensory theta was a reliable indicator of immediate tremor amplitude.
These results overall argue for a combined subcortical–cortical stimulation/recording paradigm not unlike cortical–thalamic closed-loop DBS for ET (Opri et al., 2020). By combining recordings from the STN and sensorimotor cortex, an algorithm could infer whether tremor was about to emerge (STN and MC theta) or was already present (SC theta). In particular, somatosensory cortical recordings could allow for continuous monitoring of tremor despite any stimulus artifact or competing oscillations in the STN. Ideally, DBS for a patient with a mixed-motor phenotype could be optimized between STN beta for bradykinesia symptoms and SC theta oscillations for tremor.
Limitations and conclusions
Because all tremor data were quantified from patients as they were moving their upper limb during our tracking task, our tremor conditions do not reflect a pure “rest” tremor. However, as parkinsonian tremor can often emerge as patients maintain a posture or perform a task, our approach still captured meaningful aspects of tremor. Because of our PD population receiving mostly STN DBS for clinical reasons, we were unable to assess the role of the GPi and cerebellar thalamus (VIM) neurophysiology to tremor onset and/or maintenance. In addition, as increased cognitive load has been found to exacerbate tremor, our observed tremor-related changes in nontremor frequencies within cortex may have reflected cognitive or visuomotor processes (e.g., eye movements) not directly related to tremor (Dirkx et al., 2020). Although we attempted to overcome the influence of individual subjects in our tremor epoch datasets by using linear mixed models, we were unable to apply linear mixed models to our directed connectivity analyses (PSI, GP) and thus may be susceptible to individual subject influence. However, our directed connectivity results were often reinforced by nondirected measures of functional connectivity (PLV), suggesting that directed results reflected the same underlying phenomena. Nevertheless, our awake behaving intraoperative recordings revealed that the STN and motor cortex work together to initiate tremor, and tremor is in part sustained by premotor–parietal synchrony.
Footnotes
This work was supported by a National Institutes of Health (NIH) National Institute of Neurological Disorders and Stroke Training Grant (T32-MH-020068) to P.M.L.; a Doris Duke Clinical Scientist Development Award (#2014101) to W.F.A.; an NIH COBRE Award from the National Institute of General Medical Sciences Grant P20-GM-103645 (Principal Investigator, Jerome Sanes) supporting W.F.A.; a Neurosurgery Research and Education Foundation grant to W.F.A.; the Lifespan Norman Prince Neurosciences Institute; and the Brown University Robert J. and Nancy D. Carney Institute for Brain Science. Part of this research was conducted using computational resources and services at the Center for Computation and Visualization at Brown University, with funding provided by NIH Office of the Director Grant S10-OD-025181. W.F.A. has received proprietary equipment and technical support for unrelated research through the Medtronic external research program. We thank the patients in this study for participation in the study. We also thank Kelsea Laubenstein-Parker for technical assistance, Karina Bertsch for administrative support, and Ann Duggan-Winkle for clinical support. In addition, we also thank Minkyu Ahn, David Segar, Tina Sankhla, and Daniel Shiebler for helping in the development of the motor task experiment.
The authors have patents and patent applications broadly relevant to Parkinson's disease (but not directly based upon this work). W.F.A. has received proprietary equipment and technical support for unrelated research through the Medtronic external research program. The authors declare no other competing financial interests.
- Correspondence should be addressed to Peter M. Lauro at peter_lauro{at}brown.edu or Wael F. Asaad at wfasaad{at}alum.mit.edu
