Cross-frequency coupling between gamma oscillations and deep brain stimulation frequency in Parkinson’s disease

Abstract

The disruption of pathologically enhanced beta oscillations is considered one of the key mechanisms mediating the clinical effects of deep brain stimulation on motor symptoms in Parkinson’s disease. However, a specific modulation of other distinct physiological or pathological oscillatory activities could also play an important role in symptom control and motor function recovery during deep brain stimulation. Finely tuned gamma oscillations have been suggested to be prokinetic in nature, facilitating the preferential processing of physiological neural activity. In this study, we postulate that clinically effective high-frequency stimulation of the subthalamic nucleus imposes cross-frequency interactions with gamma oscillations in a cortico-subcortical network of interconnected regions and normalizes the balance between beta and gamma oscillations. To this end we acquired resting state high-density (256 channels) EEG from 31 patients with Parkinson’s disease who underwent deep brain stimulation to compare spectral power and power-to-power cross-frequency coupling using a beamformer algorithm for coherent sources. To show that modulations exclusively relate to stimulation frequencies that alleviate motor symptoms, two clinically ineffective frequencies were tested as control conditions. We observed a robust reduction of beta and increase of gamma power, attested in the regions of a cortical (motor cortex, supplementary motor area, premotor cortex) and subcortical network (subthalamic nucleus and cerebellum). Additionally, we found a clear cross-frequency coupling of narrowband gamma frequencies to the stimulation frequency in all of these nodes, which negatively correlated with motor impairment. No such dynamics were revealed within the control posterior parietal cortex region. Furthermore, deep brain stimulation at clinically ineffective frequencies did not alter the source power spectra or cross-frequency coupling in any region. These findings demonstrate that clinically effective deep brain stimulation of the subthalamic nucleus differentially modifies different oscillatory activities in a widespread network of cortical and subcortical regions. Particularly the cross-frequency interactions between finely tuned gamma oscillations and the stimulation frequency may suggest an entrainment mechanism that could promote dynamic neural processing underlying motor symptom alleviation.

Video Abstract

gamma oscillations, cross-frequency coupling, deep brain stimulation, source analysis, volume of tissue activated

Introduction

Deep brain stimulation (DBS) is a well-established treatment for alleviating motor symptoms in Parkinson’s disease (Deuschl et al., 2006; Allert et al., 2018). Parkinson’s disease affects the function of cortico-basal ganglia brain networks and it has been hypothesized that DBS not only acts locally in the targeted nucleus, but that the disruption of pathological signals passing through the latter might allow physiological activity to reappear within the network’s subregions (Helmich et al., 2013; Chiken and Nambu, 2016; Muthuraman et al., 2018a).

For example, beta band oscillations, as observed in EEG, magnetoencephalography (MEG), and local field potential (LFP) recordings are increased in Parkinson’s disease patients, which might constrain neural activity into an inflexible pattern that limits information coding capacity and prevents dynamic processing within the cortico-basal ganglia network (Brown et al., 2001; Litvak et al., 2012; Brittain and Brown, 2014; de Hemptinne et al., 2015; Swann et al., 2016; Tinkhauser et al., 2017). Dopaminergic medication as well as DBS modifies beta oscillations within the basal ganglia and cortex; the extent of the observed alterations correlates with clinical improvement (Kühn et al., 2006; Neumann et al., 2016; Oswal et al., 2016; Tinkhauser et al., 2018) as well as the performance of a pronation-supination motor task (Kühn et al., 2008). Recently it was suggested that adaptive DBS can improve clinical efficiency by selectively targeting prolonged beta bursts instead of reducing overall beta power (Tinkhauser et al., 2017).

Similarly, gamma oscillations (i.e. oscillations of frequencies >30 Hz) in the basal ganglia and cortex are modulated by DBS and dopaminergic treatment. Contrary to the antikinetic effect of beta oscillations, gamma activity supposedly acts in a prokinetic fashion (Brown, 2003; Litvak et al., 2012). In the subthalamic nucleus (STN), finely tuned gamma (FTG) activity between 60 and 90 Hz is increased following dopaminergic treatment (Alonso-Frech et al., 2006; Androulidakis et al., 2007) and the extent of amplification correlates negatively with motor impairment (Lofredi et al., 2018). Litvak et al. (2012) showed that dopaminergic treatment increases movement-related gamma reactivity within the motor cortex (M1) and STN and that the extent of reactivity within the STN as well as the coherence between STN and M1 correlate with motor improvement. Gamma information flow between cortical motor regions and the STN is also increased by dopaminergic treatment (Lalo et al., 2008). However, a recent study demonstrated that STN-DBS increases gamma activity over frontal and parietal areas in Parkinson’s disease patients, which correlates negatively with symptom alleviation (Cao et al., 2017). Exaggerated cortical narrowband gamma was associated with dyskinesia and could serve as a control signal for adaptive DBS (Swann et al., 2016, 2018).

Given the proposed antagonistic functions of beta and gamma activities we raise the question whether a reciprocal linear relationship between the two frequency bands exists. Furthermore, while it was shown that adaptive DBS can act to truncate exaggerated beta bursts (Tinkhauser et al., 2017), we investigate whether modifications of gamma band activity are secondary to reduced beta activity (i.e. if reductions of beta activity allow gamma oscillations to reappear in the network) or directly relate to the stimulation itself, possibly due to an immediate entrainment (Buzsaki and Wang, 2012).

Recently, it has been proposed to study Parkinson’s disease on a systems or network-level, considering the anatomical and functional interactions of the basal ganglia with the cortex, cerebellum (CER), and brainstem (Helmich et al., 2013; Caligiore et al., 2016). Since effective DBS simultaneously affects multiple regions that are connected with the stimulation site (Koirala et al., 2016, 2018; Muthuraman et al., 2017), our understanding is that we must study the effects of DBS on a wide network instead of (pre-) selected regions. Recent advances in EEG source reconstruction have made it possible to non-invasively reveal cortical and subcortical sources with an improved spatial resolution, adding to the advantage of having a good temporal resolution (Litvak et al., 2011; Muthuraman et al., 2012, 2018b; Tamas et al., 2018; Seeber et al., 2019). We used high-density (HD)-EEG with 256 channels to investigate DBS-dependent modulations of beta and gamma oscillations during rest. The relevant regions were defined with a coherent source beamforming algorithm (Muthuraman et al., 2018b) that used the volume of tissue activated (VTA) by DBS (Horn and Kühn, 2015) as a seed region. We hypothesized that clinically effective high frequency DBS of the STN modifies beta and gamma oscillations in a wide cortico-subcortical network of connected regions and investigated whether these two bands are related by performing correlation analyses between beta- and gamma band power. To investigate whether such dynamics are coupled to the occurrence of stimulation we further calculated the cross-frequency coupling (CFC) between the stimulation frequency and beta- or gamma-band frequencies.

Materials and methods

Patient demographics and data acquisition

Thirty-one patients with clinically definite Parkinson’s disease, as diagnosed by the London Brain Bank criteria (Hughes et al., 1992), participated in our study. All patients were chronically treated with DBS for 6–12 months prior to study participation. Over this period, subjects were maintained at stimulation frequencies that were individually optimized and then fixed, irrespective of their future potential participation in the study. Sixteen patients achieved optimal improvement of motor symptoms with 160 Hz DBS. The other 15 patients received 130 Hz stimulation for optimal benefit. Clinical parameters were obtained by experienced clinicians who were blinded to the applied DBS frequency. Patient demographics are shown in Table 1. We recorded 10 min of resting state EEG activity for four different conditions: (i) stimulation switched off (Off); (ii) stimulation at clinically effective frequency (On-clinical; could be 130 Hz or 160 Hz); (iii) stimulation 20 Hz below the clinically effective stimulation (On-low); and (iv) stimulation 20 Hz above the clinically effective stimulation (On-high). Stimulation was delivered in a monopolar fashion. The field isolation containment system with front-end filters that usually compensate for high frequency magnetic resonance gradient artefacts prevented the amplifier of the MRI compatible EEG system from saturating. The order of the four conditions was pseudorandomized and we allowed for 45 min of washout period after switching the stimulation off. All participants were seated in a comfortable chair in a slightly reclined position, with both forearms supported to the wrists by firm armrests. Dopaminergic medications were withdrawn at least 12 h before the study.

Table 1

Demographics and clinical details of study participants

Parameters	Stimulation 160 Hz	Stimulation 130 Hz
n	16	15
Male/female	11/5	9/6
Mean age, years	65.4 ± 4.8	68.8 ± 4.7
Disease duration, years	15.5 ± 2.4	14.46 ± 2.9
LED, mg
Pre	875.8 ± 175.9	932.6 ± 241.5
Post	351.4 ± 111.3	355.3 ± 117.6
UPDRS III MED OFF DBS Off	35.8 ± 3.2	35.9 ± 4.6
UPDRS III MED ON DBS On	20.5 ± 3.4	22.4 ± 4.0
UPDRS III MED OFF DBS On-low	35.7 ± 3.3	37.2 ± 4.2
UPDRS III MED OFF DBS On-high	35.6 ± 3.4	34.6 ± 3.2

Details are provided separately for the two groups based on the stimulation frequency. LED = levodopa equivalent dosage; MED OFF = without medication; UPDRS = Unified Parkinson's Disease Rating Scale.

EEG was recorded with a high density 256-channel recording system [Electrical Geodesics, Inc. (EGI), Philips], with Cz as reference at a sampling rate of 5000 Hz. Data were analysed offline. All patients underwent preoperative MRI using a 3 T MRI scanner (Siemens TrioTim) with a 32-channel head coil and postoperative CT scan using a Toshiba Aquilian with slice thickness of 0.5 mm. This included whole-brain high resolution T₁-images using standard MPRAGE (magnetization-prepared 180° radio-frequency pulses and rapid gradient-echo) sequence with repetition time = 1900 ms, echo time = 2.52 ms, flip angle = 9 and voxel resolution of 1 × 1 × 1 mm³. The study protocol was approved by the local ethics committee in Mainz and all patients provided written informed consent before the procedure.

Preprocessing and time frequency analyses

EEG data preprocessing and part of the spatial filter analyses were performed using MATLAB (version 2015a, Mathworks Inc.) and the fieldtrip toolbox (Oostenveld et al., 2011). The initial preprocessing steps were performed by a researcher who was blinded to the stimulation conditions. Initially, EEG data were re-referenced to the common grand average reference of all channels. The raw data were low-pass filtered (fourth-order Butterworth filter; cut-off frequency: 500 Hz) to avoid aliasing, followed by high-pass filtering at 0.5 Hz. Then, data were subjected to independent component analyses (FastICA) to remove components related to DBS, muscle, eye blink, eye movement and line noise artefacts. On average, 16 of 256 components [16 ± 2.3, mean ± standard deviation (SD) were rejected (DBS artefact: 6 ± 1.24; eye artefact: 5 ± 0.68; line noise: 2 ± 0.34; muscle artefacts: 2 ± 1.21]. Residual muscle artefacts were visually inspected, removed and interpolated with the cubic interpolation method. Continuous data were then decomposed into their time-frequency representation using the multitaper method (Mitra and Pesaran, 1999; Muthuraman et al., 2010a). Seven orthogonal tapers were used with good leakage and spectral properties, and discrete prolate spheroidal sequences (DPSS) were applied (Pollak and Slepian, 1961).

Reconstruction of brain activity

To solve the EEG forward problem in source reconstruction, we estimated a volume conduction model and the lead-field matrix (LFM) containing information about the geometry and conductivity with the finite-element method (FEM) (Wolters et al., 2007). Skin, skull, CSF, and grey and white matter surfaces were extracted from the individual anatomical T₁-MPRAGE scans and individual electrode locations were used. A complete description has been described previously (Muthuraman et al., 2010b, 2012). A full description of the beamformer linear constrained minimum variance spatial filter is given elsewhere (Van Veen et al., 1997; Muthuraman et al., 2018b). The output of the beamformer at a voxel in the brain can be defined as a weighted sum of the output of all EEG channels. The frequency components and their linear interaction are represented as a cross-spectral density (CSD) matrix. To visualize power at a given frequency range, we used a linear transformation based on a constrained optimization problem, which acts as a spatial filter (Van Veen et al., 1997). The spatial filter assigned a specific value of power to each voxel. For a given source the beamformer weights for a location of interest are determined by the data covariance matrix and the LFM. Voxel size was 2 mm, resulting in 6676 voxels covering the entire brain.

Network definition and VTA analysis

To identify distributed regions that are affected by DBS we used a coherent source beamforming algorithm as described in previous work (Muthuraman et al., 2018b). In short, this method uses a reference signal to identify other sources displaying increased coherence in an iterative manner. This reference signal was individually extracted from a seed region within the STN that is targeted by DBS as estimated by the VTA by stimulation. This was done based on the clinically optimal DBS settings using Lead-DBS [a MATLAB based toolbox for reconstructing the implanted electrodes and simulating the stimulations (https://www.lead-dbs.org/)]. The details of electrode localization and VTA reconstruction have been previously described in detail (Horn and Kühn, 2015; Horn et al., 2019). In brief, preprocessing was performed using SPM12 (http://www.fil.ion.ucl.ac.uk/spm/software/spm12) where postoperative images were linearly co-registered to preoperative MRI and were manually controlled for each patient and refined if needed. Obtained images were then normalized into ICBM 2009b NLIN asymmetric space based on the preoperative MRI and finally DBS electrode contacts were localized within MNI space using Lead‐DBS software. To construct a model for volume conduction of the DBS electrode from the active contact, a tetrahedral volume mesh was generated based on the surface meshes of DBS electrodes and subcortical nuclei using the Iso2Mesh toolbox (http://iso2mesh.sourceforge.net/) as included within Lead‐DBS. All parameters used for the reconstruction were as previously published in Horn et al. (2019) and the voltage applied to the active electrode contacts was introduced as a boundary condition (Åström et al., 2009).

Individual masks were created by taking the pooled time series from the individual VTA voxels as a reference signal. For the beta (14–30 Hz) and gamma bands (31–100 Hz), a within-subject surrogate analysis was used to define the significance level to identify activated voxels in other regions. Their activity was then extracted from the source space. In a further analysis, all the original source signals for each region with several activated voxels were combined by estimating the second order spectra and employing a weighting scheme depending on the analysed frequency range to form a pooled source signal estimate for each region as previously described (Rosenberg et al., 1989; Amjad et al., 1997; Muthuraman et al., 2014). This network was found iteratively starting from the strongest coherent source to the weakest source. We found the sources in three cortical regions [M1, premotor cortex (PMC), supplementary motor area (SMA)] and two subcortical regions (STN, CER) for each patient separately. Additionally, the posterior parietal cortex (PPC) was used as a control region for all our analyses; the MNI coordinates [–41, –67, 40] were taken from our previous work on parkinsonian patients (Muthuraman et al., 2018b). The regions’ individual subject MNI coordinates were then also used for the On-clinical condition and the other three conditions (Off, On-low and On-high) as reference.

Spectral analysis and cross-frequency coupling

We separately analysed the power over five different frequency bands, namely delta (1–3 Hz), theta (4–7 Hz), alpha (8–13 Hz), beta (14–30 Hz) and gamma (31–100 Hz) in all six regions. The absolute power was analysed from both ipsilateral and contralateral hemisphere and pooled together for all four conditions separately.

The cross-frequency inter-regional coupling measure of power-power indicates how amplitude modulations in one frequency in one region depend on amplitude modulations in another frequency and region. In contrast to coherence measures, the advantage of this method is that couplings can be detected between different frequencies and even within the same region. The estimation is done by comparing the envelopes of the signals at the given frequencies from the two regions. We estimated the relationship between the amplitude of the higher frequency signals (including the stimulation frequency) and the low frequency signals (i.e. 1–100 Hz) by correlation. The CFC was estimated with a window length of 5 s with a 50% overlap. The complete analysis pipeline is shown in Fig. 1.

Figure 1

Schematic of the analysis pipeline. (A) In a first step the VTA was calculated using lead-DBS (see also Fig. 2). (B) Utilizing the reconstructed source time series of voxels identified by VTA analysis as a reference for the clinical stimulation condition we identified six coherent sources that were used for further analyses [STN, CER, M1, premotor cortex (PMC), and supplementary motor area (SMA)]. Additionally, the PPC was included as a control region. PD = Parkinson’s disease.

Statistical analyses

Differences in pre- and post-surgery medication were evaluated by paired t-tests. To investigate the effect of the four stimulation conditions on motor impairment, a one-way repeated measures ANOVA with the within subject factor ‘stimulation condition’ (Off, On-clinical, On-low, On-high) was performed with Unified Parkinson’s Disease Rating Scale – Part 3 (UPDRS III) scores as dependent variable. Similar ANOVA’s were performed for the source power of different frequency bands within the identified regions to investigate the effect of the four stimulation conditions on oscillatory activity within the network. Post hoc tests for all ANOVA’s were Bonferroni corrected.

To examine whether the amount of gamma power modulation is related to the amount of beta power modulation, we performed Pearson correlation analyses between beta and gamma power for the Off and On-clinical conditions and the difference between On-clinical and Off (On-clinical − Off) of beta (dBeta) and gamma power (dGamma).

The significance of power-to-power CFC was determined by surrogate analyses based on the bootstrapping method (Kamiński et al., 2001). Time series were shuffled 1000 times and frequency-frequency clusters showing enhanced CFC in the original time series compared to random shuffling were identified. Additionally, one-way repeated measures ANOVA with the within subject factor ‘epoch’ of 1 min duration was performed separately for the three conditions (On-clinical, On-low, On-high) to investigate the temporal dynamics of CFC (power-to-power) over M1. Pearson correlation was then performed to investigate the association between CFC and the UPDRS III scores in the On-clinical condition.

Data availability

The data that support the findings of this study are available from the corresponding author, upon reasonable request.

Results

Motor improvement by clinically effective DBS

To confirm that only effective DBS alleviates motor impairment, UPDRS III scores were compared for all four conditions. Mean scores were 35.90 ± 3.90 (Off), 21.48 ± 3.80 (On-clinical), 36.48 ± 3.98 (On-low), and 35.12 ± 3.33 (On-high). One-way repeated measures ANOVA showed a statistically significant effect of stimulation condition on UPDRS III scores [F(3,120) = 116.30, P < 0.0001, ηp2 = 0.25]. Post hoc analysis revealed that only clinically effective stimulation (On-clinical) significantly (P < 0.0001) reduced UPDRS III scores. The l-DOPA equivalent dosage was significantly decreased 6 months after the DBS procedure (353.32 ± 112.56 mg, post-DBS) in comparison to pre-DBS (903.35 ± 208.68 mg, pre-DBS) (t = 12.12, df = 30, P < 0.0001).

VTA results

We did not find a significant difference (t = 0.86; P = 0.3276) between the volumes of the VTA in the high (160 Hz) and low frequency (130 Hz) groups. VTA-STN overlap was not significantly different between the low and high frequency group based on the Euclidean distance from the centre of VTA MNI coordinates (t = 0.94, P = 0.2435), demonstrating that DBS stimulates similar tissue in both groups. The individual electrodes and clinical contact based on the optimal clinical frequency are shown in Fig. 2A and B, respectively. The VTA for all the subjects are shown in Fig. 2C. No significant correlation was found between the volumes of VTA and age, sex, or disease duration.

Figure 2

Electrode location and VTA. (A) DBS electrode reconstruction for all subjects using lead-DBS. Red and white colours represent electrodes of patients stimulated with 130 Hz and 160 Hz, respectively. Subcortical structures are based on DISTAL atlas [orange: STN, green: internal globus pallidus (GPi), blue: external globus pallidus (GPe), and red: red nucleus] laid over a 7 T MRI ex vivo 100-μm thick human brain background template. Electrode contacts used for stimulation are shown in green (130 Hz) and magenta (160 Hz) and (B) VTA is shown in cyan (130 Hz) and magenta (160 Hz). (C) All electrodes locate to the sensorimotor region of the STN and VTA simulation shows that DBS primarily activates tissue in this area.

Coherent cortical and subcortical sources

Using the individual VTA pooled time series as the reference, we identified coherent sources in the regions of M1, PMC, SMA, STN and CER. All of the identified sources were statistically significant for the beta band (M1: t = 12.46, P < 0.0001; PMC: t = 9.64, P < 0.0001; SMA: t = 12.42, P ≤ 0.0001; STN: t = 8.98, P < 0.0001; CER: t = 9.93, P < 0.0001) and for the gamma band (M1: t = 11.23, P < 0.0001; PMC: t = 10.25, P < 0.0001; SMA: t = 11.12, P ≤ 0.0001; STN: t = 9.24, P < 0.0001; CER: t = 8.42, P < 0.0001). These regions were then examined in the further analyses. The PPC was added as a control region.

Effective DBS reduces beta power and increases gamma power

To identify the effect of DBS on resting state oscillatory activity, one-way repeated measures ANOVAs were performed for five different frequency bands (delta, theta, alpha, beta, and gamma) in each of the identified regions separately to compare source power for all four stimulation conditions (Off, On-clinical, On-low, On-high). The results for the beta and gamma bands are presented in Fig. 3.

Figure 3

Stimulation induced beta- and gamma-band power changes. One-way repeated measures ANOVAs revealed that clinical stimulation significantly reduced beta-band and concurrently increased gamma-band power in M1, PMC, SMA, STN, CER, but not in PPC (Post hoc test with Bonferroni procedure, all P < 0.05). Stimulating with higher or lower frequencies (in comparison to clinical frequencies) exclusively reduced beta power in the STN.

There were statistically significant main effects of stimulation on spectral power in the beta and gamma bands in M1 [F(3,120) = 279.53, P < 0.0001, ηp2 = 0.25; F(3,120) = 81.14, P < 0.0001, ηp2 = 0.25], PMC [F(3,120) = 90.94, P < 0.0001, ηp2 = 0.25; F(3,120) = 17.96, P < 0.0001, ηp2 = 0.25], SMA [F(3,120) = 90.94, P < 0.0001, ηp2 = 0.25; F(3,120) = 17.96, P < 0.0001, ηp2 = 0.25], STN [F(3,120) = 79.23, P < 0.0001, ηp2 = 0.25; F(3,120) = 255.16, P < 0.0001, ηp2 = 0.25], and CER [F(3,120) = 156.76, P < 0.0001, ηp2 = 0.25; F(3,120) = 22.55, P < 0.0001, ηp2 = 0.25]. Post hoc analyses (Bonferroni procedure, all P < 0.05) revealed that only clinical stimulation (On-clinical) significantly reduced beta power and increased gamma power. Interestingly, however, beta power was also reduced for On-low and On-high in the STN. Spectral power in the delta, theta, and alpha bands was not significantly affected by any condition in any region. Additionally, there was no statistically significant effect of stimulation on spectral power in any of the bands in the PPC [delta: F(3,120) = 1.52, P > 0.05; theta: F(3,120) = 1.91, P > 0.05; alpha: F(3,120) = 0.24, P > 0.05; beta: F(3,120) = 0.60, P > 0.05; gamma: F(3,120) = 0.61, P > 0.05]. The results of these tests are summarized in Supplementary Table 1. To summarize, ANOVAs revealed that only the clinically effective stimulation could reduce beta power and increase gamma power in M1, PFC, SMA, and CER, but not in the PPC. The STN was the only region that showed reduced beta power in all three stimulation conditions, while gamma power was again exclusively increased for the On-clinical condition.

In the On-clinical condition, there was a significant negative linear correlation between beta- and gamma power for all regions except the PPC (M1: r = −0.5616, P = 0.001; PMC: r = −0.5147, P = 0.0045; SMA: r = −0.5078, P = 0.0064; STN: r = −0.5279, P = 0.0023; CER: r = −0.4984, P = 0.0048, PPC: r = 0.1472, P = 0.724). Similar correlations were found between dBeta and dGamma (M1: r = −0.4470, P = 0.017; PMC: r = −0.4754, P = 0.0057; SMA: r = −0.4587, P = 0.0078; STN: r = −0.4189, P = 0.0043; CER: r = −0.4374, P = 0.024, PPC: 0.1026, P = 0.845). When the stimulation was switched off there was no significant correlation between beta and gamma band power. Figure 4 shows correlations for M1 and STN for both absolute power in the On-clinical condition and the difference between the Off and the On-clinical conditions. All P-values are Bonferroni-corrected.

Figure 4

Association between beta and gamma power. Correlation analysis demonstrated a significant negative correlation of beta and gamma power during the clinically effective stimulation condition (On-clinical) for M1 and STN. In addition, the difference between the On-clinical and Off condition of beta and gamma power was significantly negatively correlated. The coloured regions denote the 95% confidence intervals for the correlation. dBeta and dGamma refer to (On-clinical − Off) power estimates. The r² and P-values are included separately for each subplot.

Narrowband gamma activity is coupled to VTA effects at stimulation frequency

To investigate whether any of the DBS-related changes in oscillatory activities were directly coupled with stimulation driven changes in the VTA region, we calculated the power-to-power CFC between frequencies from 1 Hz to 100 Hz in the regions identified above and the power between frequencies from 1 Hz to 200 Hz of the signal derived from the VTA (Fig. 5). In all regions that exhibited beta and gamma modifications, statistically significant clusters were observed between a narrow gamma band with frequencies between 60 Hz and 80 Hz and the power at the stimulation frequency in the VTA (Supplementary Table 2). We did not find any significant clusters for the On-low (110 and 140 Hz) and On-high (150 and 180 Hz) conditions, as shown in Supplementary Figs 1–4. We found significant temporal change in the CFC (power-to-power) in M1 for the factor ‘epoch’ [F(54,2987.5) = 6.74, P < 0.0001, ηp2 = 0.108] only for the clinically effective frequency (130 and 160 Hz) (P < 0.0001) as shown in Fig. 6A. Epoch 1 was significantly different to all the later epochs (P < 0.0001). Epochs were not significantly different for the On-low and On-high conditions (P > 0.05) as shown in Fig. 6B and C. This demonstrates that the observed CFC under the On-clinical condition needs ∼1 min to establish. Importantly, such delay illustrates that the CFC does not relate to a stimulation artefact, since this would be established with the onset of stimulation. Similar dynamics have also been noted with another stimulation-related phenomenon, the event-related resonance activity recorded at the stimulation site in the STN. Here the frequency and amplitude of the event-related resonance activity also takes ∼1 min to reach steady state (Wiest et al., 2020). Furthermore, there was a significant negative correlation of the gamma CFC in M1 with the UPDRS III in the On-clinical condition (r = −0.5266; P = 0.0023) as shown in Fig. 7.

Figure 5

Inter-regional CFC between the VTA reference signal and the regions of interest. CFC indices revealed clusters of narrowband gamma power coupled to the stimulation frequency (A) 130 Hz and (B) 160 Hz, respectively. The black line in the colour bars indicates the significance threshold based on surrogate analyses. See Supplementary Table 2 for ANOVA results and exact P-values.

Figure 6

Dynamics of CFC in M1. The temporal variation of CFC (power-to-power) in M1 for the entire 10 min divided into 10 epochs of 60 s. Stimulation was switched on at time point 0. The CFC is only significantly increased during On-clinical stimulation and takes ∼1 min to establish (A). This demonstrates that increased CFC does not relate to a stimulation artefact. If this was the case, increased CFC would be visible in non-clinical stimulation and would develop from the onset of the stimulation. (A) On-clinical 130 Hz (orange) and 160 Hz (blue); (B) On-low 110 (orange) and 140 Hz (blue); and (C) On-high 150 Hz (orange) and 180 Hz (blue). Statistical analyses between the epochs are shown with red lines: epoch 1 is significantly different to epochs 2 to 10, as shown by the black lines, but only for the On-clinical frequencies 130 Hz and 160 Hz (P < 0.0001). The bold line indicates the mean and the shaded region indicates the standard deviation over all subjects. Note the y-axis in A differs from that in B and C.

Figure 7

Association between gamma CFC and UPDRS III. Correlation analysis demonstrated a negative correlation of gamma CFC and UPDRS III On-clinical scores during clinically effective stimulation (r = −0.5266; P = 0.0023).

Discussion

This study shows the power-to-power CFC of narrow band gamma oscillatory activity within multiple regions of the cortico-basal ganglia network to power at the frequency of clinically effective DBS recorded in the VTA during rest. We further demonstrate a negative correlation between beta and gamma band power, and between the extent of cross frequency coupling of narrow-band gamma with effective stimulation frequency and clinical impairment, highlighting the clinical relevance of DBS induced power-to-power-coupling. These findings suggest that entrainment of FTG oscillations might promote dynamic neural processing, facilitating improvements of motor function in patients with Parkinson’s disease.

DBS alleviates motor symptoms and modulates narrow band gamma activity

We observed a power-to-power CFC between narrowband gamma oscillations in regions far from the stimulation target and stimulation frequency activities in the VTA. This was only true for stimulation at clinically effective DBS, was not seen in the control posterior parietal region and was delayed in onset, taking ∼1 min to establish itself. These features collectively make it unlikely that the coupling was related to local or volume conducted stimulation artefacts. Moreover, the narrow band gamma power-to-power CFC itself correlated negatively with motor impairment. Frequencies >30 Hz may relate to several independent processes and it is important to distinguish between broadband (30–200 Hz) and FTG activity within a narrow band between 60–90 Hz, thought to be responsible for arousal, attention, or vigour-related functions (Litvak et al., 2012; Jenkinson and Kuhn, 2013). Fischer et al. (2017) demonstrated that successful stops in a modified stop-signal task were preceded by narrowband gamma activities before the patients knew they had to stop, suggesting that FTG relates to facilitation of dynamic local processing rather than being uniformly prokinetic. In this sense the FTG could be viewed as the opposite of beta activity, which promotes the status quo (Gilbertson et al., 2005; Engel and Fries, 2010), and it is interesting to note the reciprocal relationship between these activities in several key motoric regions in this study. Interestingly, stimulating the STN with individually detected movement and medication-dependent gamma peak frequencies within the range of FTG frequencies produced similar clinical benefits as stimulating with frequencies >100 Hz (Tsang et al., 2012). In addition, FTG may spontaneously involve harmonic components that lie in the range of clinically effective DBS frequencies, in the absence of DBS (see Fig. 2 in Brown et al., 2001).

The observed power-to-power CFC between narrow band gamma oscillations and the stimulation frequency VTA effects in our study suggest that DBS with clinically effective frequencies may trigger intrinsic FTG oscillations via an entrainment mechanism (Buzsaki and Wang, 2012; Li et al., 2012; Agnesi et al., 2015; Swann et al., 2016), potentially modifying neuronal communication in a temporally coordinated manner, as proposed by Fries (2015) and thereby alleviating motor symptoms in Parkinson’s disease. Swann et al. (2016) have previously demonstrated that the florid FTG recorded at the cortical level ON medication may be entrained by high frequency stimulation of the STN when patients experience dyskinesia. This entrainment was manifest as a shift in the frequency of the FTG peak to that of the subharmonic (65 Hz) of the stimulation frequency (130 Hz). In contrast, CFC did not occur at precisely subharmonic frequencies in our data, instead generally occurring at a frequency that fell between the subharmonics of the two DBS frequencies used. However, we applied stimulation when patients were OFF levodopa, when endogenous gamma oscillations are known to be greatly attenuated and less well tuned (Jenkinson and Kuhn, 2013). These characteristics suggest that gamma oscillations in the OFF medication state may be over-damped and therefore less susceptible to frequency shifting when forced by external stimulation in the form of DBS. The present findings raise the possibility that OFF medication DBS is also able to entrain FTG, and that these effects may be mediated by changes in the VTA/STN region driven by clinically effective stimulation. One candidate phenomenon in the STN is the event-related resonance potential (ERNA), which is harmonically related to the stimulation frequency and takes ∼1 min to develop (Sinclair et al., 2018, 2019).

The fact that the observed power-to-power CFC was evident both in the absence of a specific task or dyskinesia and correlated negatively with motor impairment suggests that entrainment of such activities might promote dynamic processing, in line with the proposed function of FTG oscillations (Litvak et al., 2012). The absence of power-to-power CFC between beta band activity and the stimulation frequency, on the other hand, indicates that beta power is not linearly related to any direct effects of conventional DBS.

DBS modulates cortical and subcortical beta and gamma activities

Clinically effective DBS simultaneously reduced beta band power and amplified gamma band power in all the regions that were identified by coherent source beamforming of the VTA time series. Interestingly, although beta power in the STN was also reduced by On-low and On-high stimulation, this was not the case in any other region, suggesting that motor symptom alleviation cannot be achieved by exclusively reducing STN beta power but critically depends on network-wide effects. Beta-band oscillations are likely related to maintaining a current sensorimotor or cognitive state (Engel and Fries, 2010) and exaggerated beta oscillations in Parkinson’s disease patients might limit the information coding capacity, leading to inflexible processing (Brittain and Brown, 2014; Tinkhauser et al., 2018). Finding attenuated beta power within the STN is in line with recent observations that DBS alleviates excessive local synchronization in this frequency band (Kühn et al., 2008; Oswal et al., 2016; Tinkhauser et al., 2017). However, a stimulation-dependent reduction of this rhythm within cortical regions is more controversial. While some studies show little or no modifications of cortical beta power (de Hemptinne et al., 2015; Oswal et al., 2016), others show beta modulations (Whitmer et al., 2012; Abbasi et al., 2018). Whitmer et al. (2012) found that DBS attenuated cortical signals in the beta band in close vicinity to the putative origin of the SMA/M1 hyperdirect pathway, whereas activity of less specific regions of the motor cortex remained unaffected. Such spatial specificity supports the hypothesis that cortical activity might be modified by DBS antidromically. Importantly, extensive spatial smoothing in the forward solution could impede the detection of focal power changes. In comparison to Oswal et al. (2016), who used a voxel size of 5 mm for the forward solution and did not show a significant beta-decrease, a voxel size of 2 mm used in our study seems to increase the ability to detect local beta power decreases in the presence of spatial selectivity. In a more recent MEG study, Abbasi et al. (2018) used the same source reconstruction method as we did on rest recordings to demonstrate that DBS significantly suppressed alpha/low beta (8–22 Hz) oscillations over M1, SMA, and PMC.

We further identified attenuated beta activation within the CER. In one of our recent studies, we demonstrated the pathological involvement of the CER in parkinsonian tremorgenesis (Muthuraman et al., 2012, 2018b) and revealed bidirectional interactions between the CER and M1/SMA during voluntary movements (Tamas et al., 2018). In addition, evidence that di-synaptic pathways project from the STN to the CER via the pontine nuclei has been reported in primates (Bostan et al., 2010), providing a possible anatomical link between these regions that could serve as an interface for the observed power modulations in the CER.

One possible mechanism for DBS to suppress beta oscillations not only locally within the STN but also in the cortex is the antidromic activation of deep layers within the cortex, which might supress cortical beta power (Li et al., 2007; Whitmer et al., 2012). This is supported by a study in rats showing antidromic recruitment of M1 layer V neurons with optical stimulation in STN (Gradinaru et al., 2009). An alternative explanation is related to the disruption of pathologically increased information flow passing through the STN. Untreated, excessive hyperdirect cortico-STN beta-activity might re-innervate the cortex via the basal ganglia output nuclei. DBS could increase afferent GABAergic input from the external part of the globus pallidus (GPe) to the STN (Chiken and Nambu, 2016). Such inhibitory input could damp or filter cortico-STN beta activities and prevent them from being conveyed through the STN (Brittain and Brown, 2014; Chiken and Nambu, 2016).

In addition to the DBS-induced beta power reduction, we found a concurrent increase of gamma power. Increased cortical gamma synchronization has previously been related to DBS (McCracken and Kiss, 2014; Cao et al., 2017) and dopaminergic treatment (Androulidakis et al., 2007; Litvak et al., 2012). Furthermore, we found a negative linear correlation between the power in the two bands during clinically effective stimulation which was not present when stimulation was switched off. Therefore, patients with low beta power concurrently exhibit high gamma power when stimulated. Furthermore, patients that showed greater reductions in beta power between off and On-clinical conditions also exhibited larger increases in gamma, as demonstrated by the negative correlation between dBeta and dGamma. This suggests that, while gamma and beta oscillations may act independently in the absence of stimulation, DBS might elicit a negative relationship between the two bands, thereby reinstating a balance between activities promoting dynamic processing and activities promoting the status quo. A possible explanation for the observed gamma power increase in our study is that DBS-induced attenuation of beta activity allows gamma activity to re-emerge in the subcortico-cortical network. Alternatively, the potential entrainment of FTG activity during clinically effective DBS might suppress beta activity. Either way, normalizing the balance between beta and gamma activities could improve flexible on-demand information processing during motor performance.

Significance of DBS effects on the subcortico-cortical network

To our knowledge, this is the first study demonstrating that DBS, when performed with the clinically effective frequency, concurrently affects oscillatory activity in PMC and CER alongside STN, M1, and SMA. We revealed modulations in the beta and gamma band in all of these regions. Therefore, we postulate that these network-wide DBS effects might counteract network-wide pathological phenomena to alleviate motor impairment. Particularly the gamma-band modulations could be caused by an entrainment mechanism, as suggested by high CFC with stimulation frequency power. Stimulation frequencies other than the optimal frequency might lack the ability to trigger such resonance. On the other hand, the absence of power-to-power CFC between beta and stimulation frequency indicates that beta power is not linearly related to any direct effects of conventional DBS.

The exact mechanisms and pathways by which STN-DBS affects the subcortico-cortical network cannot be ascertained from the current work and merit further investigation. It should be noted that we cannot infer any directionality or relative circuit conduction delays from our analysis, which prevents us from separating the relative contributions of antidromic activation of the cortex via the hyperdirect pathway and orthodromic activation via the pallido-thalamo-cortical feedback loop.

Source analysis and VTA modelling

Applying coherent source beamforming in conjunction with VTA modelling to resting state recordings, we illustrate DBS-dependent spectral alterations of beta and gamma oscillations over a wide network of regions that are connected to the STN. We have previously demonstrated that distinct networks involved during voluntary or involuntary rhythmic movements can be revealed by applying coherence analysis on the source level with the aid of peripheral muscular electromyography (Anwar et al., 2016; Muthuraman et al., 2018b; Tamas et al., 2018). Here we demonstrate that the source activity within modelled VTA is suitable as a reference to identify regions that exhibit DBS-dependent oscillatory modifications.

Clinical relevance

Altogether our findings provide further evidence that effective DBS modifies oscillatory activity within a distinct network of regions connected to the stimulation site. Most of the studies focusing on oscillatory activity within subcortical regions of Parkinson’s disease patients were based on LFP recordings from externalized electrodes, limiting the opportunities to collect information about the long-term and distributed effects of DBS (Neumann et al., 2019). Additionally, in the limited time interval between surgery and LFP recording, potential stun effects cannot be ruled out (Chen et al., 2006; Mann et al., 2009). In previous work we demonstrated that HD-EEG recordings can be successfully utilized to obtain signals from cortical as well as subcortical regions (Muthuraman et al., 2012, 2018b; Anwar et al., 2016; Tamas et al., 2018). Here we extend these findings by revealing cortical and subcortical DBS-induced oscillatory modifications during rest recordings. Therefore, HD-EEG proves to be an effective non-invasive tool for studying pathological and treatment-related electrophysiological adaptations in patients, even months after DBS surgery.

Efforts for finding stable and reliable biomarkers for adaptive closed-loop stimulation are currently focused on LFP recordings in the STN and motor cortex (Tinkhauser et al., 2017; Swann et al., 2018). The finding that DBS effectively modulates a variety of cortical and subcortical regions substantially increases the number of potential locations and measures that could be exploited for controlling the stimulation. In addition, in clinical practice, finding the optimal stimulation settings can constitute a time-consuming task, which is currently based on subjective clinical UPDRS ratings (Neumann et al., 2019). Although the UPDRS is a widely accepted and relatively reliable tool, inter-rater reliability might not be optimal, especially when performed by less experienced individuals (Siderowf et al., 2002). Furthermore, motor improvement evaluated unblinded can be overestimated by up to ∼20% in comparison to blinded assessment (Little et al., 2013). Hence, further advances in the development of HD-EEG are needed that include the automation of source reconstruction and algorithms to probe relevant DBS-induced electrophysiological modification. This could provide additional objective markers that assist in finding the optimal stimulation parameters.

Limitations

While we demonstrated that effective DBS improves motor performance, the fact that DBS at a frequency that was 20 Hz lower or higher had no effect on the UPDRS III scores is surprising. Nevertheless, a systematic varying of stimulation frequency in small steps has shown that only a limited number of frequencies robustly improve movement amplitude in a tapping task, and that even small frequency deviations can drastically reduce efficacy (Huang et al., 2014). Similarly, Di Giulio et al. (2019) showed peaked DBS frequency tuning curves in gait performance.

In conjunction with the fact that CFC was not significantly increased by non-clinical stimulation, the lack of clinical efficacy when changing stimulation frequency by as little as 20 Hz in our study could be related to long term effects of chronic stimulation. All participants in our study were chronically stimulated for 6–12 months. Figure 6B indicates that DBS with ineffective frequencies may slightly increase CFC after 1 min, although this is not significant. This suggests that there is a very small physiological effect of inefficient stimulation, which might not be sufficient to influence the clinical state of the patients. At this point we can only speculate that with prolonged DBS, functional adaptations within the network might induce a preferential resonance or tuning towards the applied stimulation frequency, explaining the observation that only the clinically optimized chronic DBS frequency increased CFC and improved motor performance. During the initial programming shortly after DBS surgery, such tuning might not be present, leading to clinically positive responses over a range of frequencies. Further studies are needed to investigate whether tuning is modulated by chronic stimulation. Such studies should include a systematic testing of smaller frequency intervals as well as serial testing of stimulation frequencies that begins soon after surgery.

Conclusion

This study is the first to demonstrate concurrent oscillatory modifications introduced by STN-DBS over a distributed cortico-subcortical network with the help of HD-EEG, showing that STN-DBS modifies activities in multiple connected regions. Such network-wide impact might ease network-wide pathological activities. Specifically, the power-to-power CFC between FTG and VTA power at the stimulation frequency, which correlates with motor scores, provides further impetus to study the role of gamma oscillations in normal and abnormal motor control. Future studies, including effective connectivity measures, could substantially aid in disentangling the possible mechanisms that cause the observations made here. Biomarkers derived from HD-EEG have the potential to provide an objective and reliable tool to support clinicians in finding optimal stimulation parameters.

Acknowledgements

We are grateful to Cheryl Ernest for their comments and help in proofreading the manuscript. Parts of this research were conducted using the supercomputer Mogon and advisory services offered by Johannes Gutenberg University Mainz (hpc.uni-mainz.de), which is a member of the AHRP and the Gauss Alliance e. V. We gratefully acknowledge the computing time granted on Mogon.

Funding

This work was supported by the German Research Foundation (DFG; SFB-TR-128, SFB-CRC 1193) and the Boehringer Ingelheim Fonds (BIF-03).

Competing interests

The authors report no competing interests.

Supplementary material

Supplementary material is available at Brain online.

 
• CER =

  cerebellum

 
• CFC =

  cross-frequency coupling

 
• DBS =

  deep brain stimulation

 
• FTG =

  finely tuned gamma

 
• PPC =

  posterior parietal cortex

 
• STN =

  subthalamic nucleus

 
• VTA =

  volume of tissue activated

Author notes

Muthuraman Muthuraman and Manuel Bange contributed equally to this work.