Neuronal oscillations predict deep brain stimulation outcome in Parkinson's disease

Background: Neuronal oscillations are linked to symptoms of Parkinson's disease. This relation can be exploited for optimizing deep brain stimulation (DBS), e.g. by informing a device or human about the optimal location, time and intensity of stimulation. Whether oscillations predict individual DBS outcome is not clear so far. Objective: To predict motor symptom improvement from subthalamic power and subthalamo-cortical coherence. Methods: We applied machine learning techniques to simultaneously recorded magnetoencephalogra- phy and local ﬁ eld potential data from 36 patients with Parkinson's disease. Gradient-boosted tree learning was applied in combination with feature importance analysis to generate and understand out- of-sample predictions. Results: A few features suf ﬁ ced for making accurate predictions. A model operating on ﬁ ve coherence features, for example, achieved correlations of r > 0.8 between actual and predicted outcomes. Coherence comprised more information in less features than subthalamic power, although in general their infor- mation content was comparable. Both signals predicted akinesia/rigidity reduction best. The most important local feature was subthalamic high-beta power (20 e 35 Hz). The most important connectivity features were subthalamo-parietal coherence in the very high frequency band ( > 200 Hz) and subthalamo-parietal coherence in low-gamma band (36 e 60 Hz). Successful prediction was not due to the model inferring distance to target or symptom severity from neuronal oscillations. Conclusion: This study demonstrates for the ﬁ rst time that neuronal oscillations are predictive of DBS outcome. Coherence between subthalamic and parietal oscillations are particularly informative. These results highlight the clinical relevance of inter-areal synchrony in basal ganglia-cortex loops and might facilitate further improvements of DBS in the future.


Introduction
Parkinson's disease is a common neurodegenerative disease, affecting approximately 6.1 M people worldwide [1]. Besides pharmacological agents such as levodopa, deep brain stimulation (DBS) is used for symptomatic treatment of Parkinson's disease. A common target structure for DBS is the subthalamic nucleus (STN) which is interconnected with the pallidum, the thalamus and several cortical areas via basal-ganglia cortex loops [2]. In patients with Parkinson's disease, activity in these loops is characterized by strong neuronal oscillations, synchronized across the connected structures [3].
Neuronal oscillations are closely related to Parkinsonian symptoms. STN beta oscillations (13e35 Hz), in particular, have been shown to reflect the motor state [4]. They are reduced by voluntary movement, pharmacological therapy and DBS [5e10]. High-gamma oscillations (60e90 Hz), in contrast, are a marker of dyskinesia, typically arising as a side-effect of dopaminergic therapy [11]. Tremor is associated with narrow-band oscillations at individual tremor frequency, observable throughout a distributed subcorticocortical tremor network [12,13]. Given their intricate relationship with Parkinsonian symptoms, recent studies have explored the utility of neuronal oscillations for optimizing DBS.
Most of these studies have focused on the dynamics of electrophysiological signals and used oscillations for adapting DBS to spontaneous changes of symptom severity, such as on-offfluctuations [14] or tremor [15,16]. Other studies have assessed the utility of quasi-stationary oscillatory activity for optimizing electrode placement, complementary to imaging studies on "sweet spots" in DBS [17,18]. Zaidel et al. demonstrated a positive linear correlation between DBS outcome and both the length of the oscillatory region in the dorsolateral STN and STN beta power [19].
Here, we adopted a similar, though more holistic approach for exploring the relationship between DBS outcome and neuronal oscillations. We analyzed simultaneous magnetoencephalography (MEG) and local field potential (LFP) recordings from Parkinson's disease patients with externalized leads to assess both STN oscillations and their synchrony with cortical activity. By applying machine learning techniques, we demonstrate that it is possible to predict DBS outcome for unseen patients based on their patterns of neuronal synchrony, considering many frequency bands and brain areas simultaneously.

Materials and methods
The aim of this study was to predict motor symptom reduction achieved by DBS based on band-limited STN power and STN-cortex coherence. For this purpose, we trained and evaluated a machine learning model operating on features extracted from MEG-LFP datasets, contributed by two previous studies performed at the University Hospital Düsseldorf [20,21]. Both studies recruited patients with Parkinson's disease selected for DBS of the STN according to standard clinical criteria.

Patient and measurement details
36 Parkinson's disease patients implanted with deep brain electrodes for STN DBS the day before the measurement took part with written informed consent, according to the Declaration of Helsinki. Patient details are given in Table S1 of the Supplementary Material. The experimental protocols were approved by the Ethics Committee of the Medical Faculty of Heinrich Heine University Düsseldorf (no. 3209 and 5608).
The experimental procedures have been described elsewhere [20,21]. Briefly, MEG signals were recorded by a 306-channel MEG system (Elekta Neuromag) with a sampling rate of 2 kHz (study 1) or 2.4 kHz (study 2). LFPs were recorded simultaneously using externalized leads and a mastoid reference. LFP signals were arranged into a bipolar montage offline. The cables used for externalization contained very little ferromagnetic material and did not cause major MEG artifacts. Forearm electromyograms as well as vertical and horizontal electrooculograms were recorded in addition. Patients were at rest in an upright position, with eyes open. The measurements took place after overnight withdrawal from dopaminergic medication (Med OFF). In a subset of patients, we performed additional recordings about 1 h after intake of levodopa (Med ON). Here, we analyzed the Med OFF data only.
The Unified Parkinson's Disease Rating Scale (UPDRS) part III of the Movement Disorders Society [22] was obtained by an experienced movement disorder specialist following optimization of DBS parameters. In most cases, scoring took place between 3 and 6 months after the LFP-MEG measurements (Table S1, Supplementary Material).

Data analysis
The general analysis pipeline is depicted in Fig. 1. It contained one sub-pipeline for feature extraction (Fig. 1A) and one for prediction (Fig. 1B).

Contact selection
First, we selected one electrode contact pair for each hemisphere by picking the contact used for therapeutic DBS at the time of UPDRS assessment and the closest neighboring contact in the direction of the midpoint between the most ventral and the most dorsal contact. This choice was adapted in case therapeutic DBS was bipolar or in case the initial choice included bad LFP channels, i.e. channels with strong noise/weak signal. In the former case, we selected the bipolar pair used for therapeutic DBS, and in the latter case, we took the closest neighbor of the bad channel in the direction of the electrode center. In case a group of segments was used for DBS in patients implanted with segmented leads, we first re-referenced the signal of each active segment from the original mastoid reference to the closest neighbor, computed features separately and averaged over segments. Lead localization, performed with LEAD DBS [23], confirmed correct placement for all electrodes under study ( Fig. 2A).

Data preprocessing
The data were preprocessed with the Fieldtrip toolbox [24]. LFP and MEG data underwent visual screening. Bad channels and epochs containing artifacts were discarded. The data were segmented into 2s windows (frequency resolution: 0.5 Hz) with 50% overlap.

STN power
We applied a Hanning taper and computed power for each integer frequency between 1 and 398 Hz using Welch's method. Line noise and its harmonics were eliminated by replacing values ± 2 Hz from the harmonics by surrogate values obtained by linear interpolation. The aperiodic (1/f) component was removed from the LFP power spectra using the fitting oscillations and one over f (FOOOF) algorithm [25]. This step was necessary to ensure that the predictive models operated on neuronal oscillations proper. Note that coherence, unlike power, is a normalized quantity not requiring this correction. When applying FOOOF, we adapted its parameters iteratively until a good fit was achieved, confirmed visually for every case. Since a good description of the entire spectrum was usually not achievable with a single model, we performed separate fits for the frequency ranges below 90 Hz and above 200 Hz (frequencies between 90 and 200 Hz were not analyzed here). The periodic minus the aperiodic component was retained and power was averaged within eight frequency bands of interest: delta/theta (3e7 Hz), alpha (8e12 Hz), low-beta

STN-cortex coherence
Coherence was estimated and localized once per frequency band rather than once per frequency. Using the multitaper method [26], we computed coherence at the band center frequency and applied appropriate spectral smoothing to include the entire band. For bands covering line noise harmonics, we computed estimates for sub-bands, excluding the harmonics, and averaged them. Coherence was source-localized using Dynamic Imaging of Coherent Sources [27]. We made use of realistic, single-shell head models based on the individual, T1-weighted MR image. The beamformer grid contained 567 locations spread out evenly across the cortical and cerebellar surface. It was aligned to Montreal Neurological Institute (MNI) space, allowing for grid parcellation into 30 supersets of regions defined in the Automatic Anatomic Labeling (AAL) atlas [28]. Details on these regions are provided in Table S2 of the Supplementary Material.
Following feature extraction, features were arranged into a feature matrix of size N patients x N features (Fig. 1). In this matrix, each subject was represented by one column comprising both STN power and STN-cortex coherence with ipsilateral and contralateral cortical parcels for both left and right STN. Two alternative designs were also tested, but found to have inferior performance (Fig. S2 of the Supplementary Material): one with hemispheres rather than subjects as unit of observation, and one in which hemispheres were ordered according their laterality with respect to the more affected body side.

Machine learning model
For predicting motor improvement, we employed extreme gradient boosting, as implemented in the XGBoost package for Python [29]. In this framework, the target score is predicted by a sequence of decision trees assembled tree-by-tree during training. Each new tree is trained on the error made by the group assembled so far, resulting in a stepwise refinement of the prediction. XGBoost has gained popularity by winning numerous machine learning competitions and is a commonly used tool in machine learning. It appears to be particularly well suited for electrophysiological datasets [30], which are typically small, structured and noisy. In a recent study by Merk et al. [31], XGBoost outperformed linear regression and artificial neural networks in the prediction of grip force based on STN and cortical oscillations.

Feature importance analysis
Feature importance analysis seeks to describe how much an individual feature or a subgroup of features contributed to a prediction made by a machine learning model. Here, we quantified feature importance using the Python implementation of SHapely Additive exPlanations (SHAP) [32]. SHAP values are estimates of Shapley values, a concept from cooperative game theory for a fair distribution of a payout among players. Besides having a range of desirable mathematical properties, SHAP values have an intuitive interpretation: they sum to the difference between the current and the average model output. While the concept is applicable to any machine learning model, specialized versions such as TreeSHAP have been developed, optimized for tree ensemble-based models such as XGBoost [33]. Following contact selection, STN power and STN-cortex coherence were computed from the Fourier spectrum. STN power underwent 1/f-correction and was averaged within frequency bands. STN-cortex coherence was source-localized using beamforming. Each source was assigned to one of 30 cortical parcels and source coherence was averaged within parcels and frequency bands. Band-limited STN power and STN-cortex coherence formed the hemisphere feature vector. (B) Leave-one-out regression. Left and right hemisphere feature vectors were stacked vertically to form the subject feature vector. The subject feature vectors were stacked horizontally to form the feature matrix. In each iteration through the leave-one-out cycle, one subject was set aside (test set). The remaining train set was divided into 3 folds for cross-validated hyper-parameter tuning and feature selection. The test features served as input to the regression model, which predicted UPDRS III sum score reduction.

Predicting DBS outcome
DBS outcome was quantified by the difference in UPDRS III sum score Med OFF/Stim OFF -Med OFF/Stim ON, unless specified otherwise. Predictions were computed sequentially for each subject in a leave-one-out fashion, i.e. each subject served as the test set once and was part of the train set in all other iterations. In each iteration through the leave-one-out loop, features were standardized using mean and variance of the train set. Next, we selected the most important k features according to the mean absolute SHAP values computed on the train set. The train set was then sub-divided into three folds for cross-validated hyper-parameter tuning with the Hyperopt package [34]. The optimization procedure and the chosen parameters are detailed in the Supplementary Material.
Model performance was quantified by the root mean squared error (RMSE) and Pearson's correlation coefficient between the actual and the predicted DBS outcome. We further applied a null model agnostic of electrophysiology for establishing a performance baseline. The null model generated predictions of DBS outcome by averaging the outcomes of the train set.

Statistics
Significance of correlation was assessed using the pearsonr function of the scipy.stats package (two-sided test; statistic: b; significance level: 0.05). When computing the correlation coefficient repeatedly, we applied false discovery rate (FDR) correction using the Benjamini-Hochberg procedure.

Features
The spectral and spatial characteristics of the features are illustrated in Figs. 2 and 3, respectively. STN power spectra contained peaks in the alpha, low-beta and high-beta band. Individual patients showed an additional high-gamma peak (Fig. 2B). The HFO spectrum was dominated by sHFO peaks, as described previously for the medication OFF state (Fig. 2C) [35,36]. The coherence spectra contained strong alpha peaks, which were ubiquitous but most pronounced in temporal areas ipsilateral to the STN (Fig. 2D). Medial sensorimotor and adjacent areas ipsilateral to the STN additionally showed strong beta peaks, as reported by previous studies (Fig. 3) [20,37]. The coherence spectra did not contain any consistent HFO peaks (Fig. 2E), but some subjects had more coherence in this range than others (Fig. S1 of the Supplementary Material). Finally, many coherence spectra had several narrow peaks in the delta/theta range, presumably reflecting tremor, occurring at slightly different frequencies in individual patients [12,13].

Model performance
We evaluated the performance of predictive models operating either on STN-cortex coherence (connectivity models) or STN power (local models) as a function of the number of features. In order to test whether the model predicted DBS benefit or symptom severity better, we predicted both the difference between the DBS OFF and the DBS ON score (benefit) and the DBS OFF score (symptom severity) separately. Note that all UPDRSIII scores were collected several months after surgery (Table S1).
When predicting DBS benefit, both connectivity and local models outperformed the null model, which estimated DBS outcome by averaging the outcomes of the train set ( Fig. 4A; RMSE Null : 6.74). Connectivity-based models outperformed the null model even with a single feature and generally performed better than local models (avg. RMSE conn : 5.1, avg. r conn : 0.64). Local models required at least four features to achieve better performance than the null model and a significant correlation between predicted and actual DBS outcomes ( Fig. 4B; avg, RMSE local : 6.0, avg. r local : 0.42).

Feature importance
This analysis aimed at revealing the most important features for successfully predicting clinical benefit. To assess feature importance, we summed absolute SHAP values over all models contributing to the previous analysis (Fig. 4) and within categories of interest such as frequency band ( Fig. 5A and B), brain region (Fig. 5C) and hemisphere with respect to the STN (Fig. 5D). For local models, STN high-beta power was the most important feature, followed by alpha and sHFO power (Fig. 5A). Strong high-beta power indicated a good DBS outcome. Strong beta-band coherence, in Fig. 3. Source-localized STN-cortex coherence. Coherence was normalized by the spatial mean prior to averaging over hemispheres. Whole-brain images were flipped such that the hemisphere ipsilateral the STN ended up on the right side. Normalized coherence is color-coded. Theta: 3e7 Hz; alpha: 8e12 Hz; low-beta: 13e20 Hz; high-beta: 21e35 Hz; low gamma: 36e60 Hz; high-gamma: 60e90 Hz; sHFO: 200e300 Hz; fHFO; 300e400 Hz. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) contrast, was not indicative of a good outcome (Fig. 5B). Connectivity models relied mostly on fHFO, low-gamma and theta oscillations. Connectivity with parietal areas was particularly important (Fig. 5C). As expected, coherence between STN and ipsilateral cortex was more important than coherence between STN and contralateral cortex (Fig. 5D).

Combining local and connectivity features
We investigated the potential merit of combining local and connectivity features when predicting DBS benefit. To simplify interpretation, we used a fixed feature set rather than running a data-driven selection for each individual subject as above. We chose the five local and the five connectivity features with the highest overall SHAP sum (as five was the lowest number for which both the local and the connectivity models reached good performance in the previous analysis). The selected features are listed in Fig. 6B. As depicted in Fig. 6A, the best-5 connectivity model (RMSE conn ¼ 3.54, r conn ¼ 0.84, p conn ¼ 1e-14) outperformed the best-5 local model (RMSE local ¼ 5.77, r local ¼ 0.52, p local ¼ 0.001). Adding local features to the connectivity features did not improve the connectivity model further (RMSE comb ¼ 3.61, r comb ¼ 0.84, p comb ¼ 1.69e-10).   Fig. 6B illustrates feature importance and the relation between feature values and feature importance for the combined model. Connectivity features generally had a stronger impact on the prediction than local features. The strongest influence was attributed to fHFO coherence between the right STN and ipsilateral inferior parietal cortex. Low values strongly drove the prediction towards worse predicted outcomes. Indeed, when relating this feature to DBS outcome directly, we observed an approximately logarithmic relationship, i.e. outcomes dropped steeply with decreasing coherence (log(coherence)-outcome correlation: r ¼ 0.40, p ¼ 0.02). Feature importance analysis further revealed that strong coupling between the right STN and ipsilateral parietal cortex and a weak coupling between the right STN and contralateral parietal cortex in the low-gamma band drove the prediction toward a good DBS outcome. So did strong coupling between the left STN and contralateral, inferior occipital cortex in the theta band and weak coupling between right STN and temporal cortex in the low-beta band.

Feature correlations
Whereas the automatically selected set of most important local features contained frequency bands with clear peaks and an established role in Parkinson's disease pathophysiology, such as the beta and the HFO band, the set of most important connectivity features did not. To better understand the nature of these features, we investigated their correlation with all other features (Fig. 7A).
Both fHFO and low-gamma coherence with parietal cortex were strongly correlated with a large set of other features, in particular >35 Hz connectivity with the entire hemisphere. Replacing each of the best five connectivity features by their closest correlate, however, decreased performance substantially ( Fig. 7B; RMSE best ¼ 3.61, r best ¼ 0.84, p best ¼ 1.69e-10; RMSE corr ¼ 6.47, r corr ¼ 0.29, p corr ¼ 0.09), demonstrating that the selected features are not arbitrary representatives of a highly correlated feature group. Destroying true phase relationships by shuffling the LFP signals in time abolished the ability to predict clinical benefit from STN-cortex coherence, demonstrating that connectivity models relied on phase information (Fig. S3 of the Supplementary Material).

Control analyses
Given the critical role of lead placement for DBS benefit, we asked whether the most informative features might reflect distance to target. None of the best 10 features were significantly correlated with the distance of the LFP channel to a published "sweet spot" for STN DBS [38], and distance could not predict clinical benefit on its own (see Supplementary Material). Similar observations were made for a set of potential confounders, including recording duration, electrode type (segmented vs. non-segmented), days passed since recording and signal-to-noise ratio. We conclude that the success of the electrophysiological models cannot be explained by correlation with these variables.

Predicting improvement of particular symptoms
In the previous analyses, we quantified DBS outcome by the reduction of the UDPRS III sum score, representing overall motor symptom severity. To see whether local and connectivity features relate to particular and possibly different symptoms, we computed predictions for each individual UPDRS III item ( Fig. 8A and C) and for the akinesia-rigidity, the tremor and the axial subset sum score ( Fig. 8B and D) using the best-5 local and the best-5 connectivity model described above. The local and the connectivity model had a very similar RMSE profile. Both predicted the DBS-induced improvement of akinesia/rigidity best.

Discussion
We have demonstrated that is possible to predict DBS outcome from STN power and STN-cortex coherence in Parkinson's disease patients. Our results indicate that STN-cortex coherence, in particular, is a good predictor of clinical benefit.

Relation to previous studies
Few studies have made out-of-sample predictions of DBS outcome based on electrophysiological data. We know of only three studies, all of which investigated signals from within or nearby the STN, recorded during surgery [39e41]. Here, we applied a network approach, incorporating both subthalamic oscillations and their synchrony with oscillatory activity in various cortical areas. In this respect, our approach can be considered an electrophysiological pendant of discriminative tractography. This is a technique for predicting DBS outcome from structural connectivity with the volume of tissue activated (VTA), an estimate of the spatial extent of neuromodulation around the active contact. This approach and related methods have facilitated accurate predictions of DBS  outcome in many disorders, including Parkinson's disease [42], obsessive-compulsive disorder [43,44], Tourette's Syndrome [45], and Essential Tremor [46]. Although structural and functional connectivity do not necessarily contain redundant information about DBS outcome [42], structural connectivity forms the basis of functional connectivity, suggesting that both might predict the efficacy of DBS in individual patients. In this study, we show that this is indeed the case. Predicted DBS outcomes, derived from a few functional connectivity features only, were highly correlated with the actual outcomes. The null model, which predicted DBS outcome based on the average outcome in the training set, performed worse, i.e. knowing the individual electrophysiology helps improving realistic expectations of DBS outcome.

The link between DBS outcome and oscillations
Predicting clinical benefit form neuronal oscillations might work for several reasons. First, exaggerated synchronization itself has been suggested to cause Parkinsonism [47]. If oscillations were the causal process, an estimate of how strong the process is around the stimulation contact might allow for predicting the DBS effect in the individual patient.
Alternatively, successful prediction might be explained by a consistent relationship between oscillations and contact location, which is known to be a crucial factor for DBS. Indeed, beta and HFO oscillations have a characteristic spatial distribution in the STN area [48,49]. Here, we found, however, that the distance of an LFP channel to a published sweet spot for STN DBS had little predictive power and did not correlate with the most informative oscillatory features [38]. This is likely a characteristic of our sample, which includes only the contacts with optimal clinical efficacy, which, in our case, differed little with respect to placement. Presumably, location would become a crucial piece of information if one were to include other contacts further off target, e.g. the contacts not selected for DBS. We could not investigate this here because we lacked a good characterization of the clinical effect for non-selected contacts.
It is possible that the prediction relied in part on information about the ability of DBS to modulate remote cortical areas, which may be key to the clinical effect of DBS [50]. This would explain why connectivity models performed better than local models and why they could predict DBS benefit but not symptom severity. Whether stimulation can reach a relevant anatomical connection might be relevant for the clinical effect of DBS, but is not relevant for the patient's motor state off stimulation.
STN-cortex coherence carries information on functional connectivity by definition. In addition, it has recently been demonstrated to correlate with the density of reconstructed fiber tracts connecting STN and cortical areas, i.e. it also contains information on anatomical connectivity [51]. Interestingly, the same study demonstrated that STN beta oscillations might arise as a consequence of cortical input, implying that even local oscillations can relate to STN-cortex connectivity.

Feature importance
In this paper, we applied feature analysis to identify particularly informative features. This analysis yielded STN high-beta power as the most import local feature, in line with a recent classification study in the primate model of PD [52]. The finding further aligns with a recent study reporting that high-beta oscillatory activity distinguishes the STN from neighboring structures [51]. STN-cortex beta coherence, in contrast, was less relevant for prediction. This observations tallies with studies on dopamine effects, which found no relation between the degree of symptom reduction achieved by levodopa and the degree of beta coherence reduction [8,53].
While beta coherence with motor cortex did not emerge as a very important feature, coherence in other frequency bands was strongly predictive of DBS outcome, even more so than STN highbeta power. Coherence between the STN and parietal cortex at high frequencies (low-gamma and fHFO) allowed for accurate estimates of DBS efficacy. This is somewhat surprising, given that these frequency bands did not contain coherence peaks. Despite this lack of structure and a strong correlation with high-frequency coupling (>35 Hz) to other cortical areas, these features appear to be particularly relevant for DBS outcome, as they were not replaceable without harming performance and lost their predictive potential through shuffling.
The importance of parietal features points towards a clinical relevance of parieto-STN connectivity. Yet, in light of the strong correlation between features, this hypothesis requires confirmation by independent studies. A number of PET studies have reported that STN DBS leads to metabolic changes in parietal areas, evidencing that DBS modulates parietal cortex [54,55]. Further, PD patients differ from healthy controls in their resting-state BOLD signal correlation between STN and parietal cortex, suggesting pathological relevance of this connection [56]. And lastly, both STN and parietal cortex were proposed to be part of brain network for response inhibition [57,58]. This network is believed to be lateralized to the right hemisphere [59,60], which might relate to the fact that all of the important parietal connectivity features identified here were right STN features, suggesting that left and right STN might not have the same relation to clinical improvement. The fact that good predictions required data from both hemispheres further supports this hypothesis.

STN power vs. STN-cortex coherence
One research question of this study was whether STN power and STN-cortex coherence carry different information about DBS outcome. Our results do not support this hypothesis. Although local models did achieve reasonable predictions, adding local features did not improve the best-performing connectivity models, indicating that the information carried by STN power was already contained in STN-cortex coherence e even though the frequency bands were different. Similarly, when predicting single UPDRS items both local and connectivity models showed a very similar performance profile, i.e. the accurately predicted symptom reductions and the non-accurately predicted symptom reductions were the same for both feature types. In conclusion, STN-cortex coherence seems to condense more information in a lower number of features, but power and coherence do not seem to be independent sources of information.
The same analysis revealed that both STN power and STN-cortex coherence predicted the reduction of akinesia/rigidity best. This might be due to the known link between akinesia and (beta) oscillations [61], and could additionally reflect common properties of the group under study. A certain level of akinesia and rigidity reduction was common to all patients, whereas only a subset of patients showed marked tremor and/or axial symptoms. Because learning occurred across subjects, oscillation-outcome relationships shared among patients were learned the best.

Limitations
While an analysis of feature importance can reveal new insights, it should be interpreted with caution. Identifying causal relationships is generally not possible with this approach [62]. Any feature related to DBS outcome might be so via correlation with other features or unobserved variables, i.e. a feature important to the model is not necessarily important to the brain. This aspect is particularly relevant with regard to a possible pathophysiological role of high-frequency coupling, which requires further investigation.
Finally, a sample size of 36 patients is large for a MEG-LFP study (as far as we know, this is the largest sample so far) but it is not at all large for a machine learning study. Future studies should aim at including more data, potentially by fostering data sharing.

Conclusions and outlook
It is possible to predict DBS outcome based on subthalamic oscillations and their synchrony with cortical oscillations. Future studies may investigate whether this link between electrophysiology and clinical improvement can be leveraged to improve lead placement and/or contact selection.

Declaration of interests
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.