Eighteen patients who met United Kingdom Brain Bank criteria for idiopathic PD underwent bilateral STN-DBS using an image-guided image-verified technique after selection by a multidisciplinary team of specialized movement disorders neurologists and functional neurosurgeons [Table 1].

All patients underwent neuropsychological screening with a Mini Mental Status Exam (MMSE) and a structural brain MRI to rule out dementia and significant brain atrophy, respectively. All patients’ UPDRS-III improved by at least a 30% preoperatively after administration of L-DOPA. For each patient, UPDRS-III scores were recorded postoperatively at 1-year follow-up in both the medication-ON and DBS-OFF (ON-OFF) and medication-ON and DBS-ON (ON-ON) states. The postoperative ON-OFF score was felt to also represent the state in which the initial rs-fMRI was acquired because improvement of UPDRS-III with L-DOPA postoperatively was correlated with preoperative L-DOPA response (r = 0.6527 [95%CI = 0.2211 0.8731], p = 0.008). Furthermore, use of the postoperative ON-OFF score accounts for disease progression and changes in response to L-DOPA over time. As such, the effects of STN-DBS over the 1-year follow-up were quantified as the percentage of improvement between the ON-OFF state and the ON-ON state.

All patients were scanned preoperatively in the medication-ON state with a 3T-Siemens Magnetom-Trio MR-B17 with a well-padded 32-channel receive head coil. Image acquisition in the medication-ON state allowed for increased patient comfort and minimization of excessive involuntary head motion, which helps to decrease the likelihood of spurious findings in rs-fMRI (Tahmasian et al., 2015). Specifically, patients remained on their normal medication regimen and the MRI followed their morning medication dose. Although this interval varied between individual patients, motor symptoms were optimally controlled at the time of MRI acquisition for all subjects.

Multi-parameter mapping sequences were acquired for structural imaging. In brief, the whole protocol consisted of three 3D-FLASH acquisitions performed with T1, proton density (PD), and magnetization transfer (MT) weighting; these were paired with B1 (transmit field-mapping data to correct for the effect of inhomogeneous flip angles on the T1 maps) and B0 field map acquisitions (spatial resolution 1 × 1 × 1 mm³; repetition time 24.5 ms; multiple echo times; field of view 256 × 256mm; flip angle 6° [PD], 21° [T1], 6° [MT]; matrix size 256 × 256; partitions 176; total acquisition time of 26 min).

For rs-fMRI, patients were instructed to maintain a fixed gaze at a crosshair while multiecho echo planar imaging sequences were obtained in two successive acquisitions, for a total duration of 15 min (spatial resolution 3 × 3 × 2.5 mm³; repetition time 2.7 s; echo time 30 milliseconds; field of view 192 × 192 mm; flip angle 90°; 45 axial slices [2.5 mm thickness]; matrix size, 64 × 64; and a total of 512 scans).

Standard preprocessing of structural and functional MRI volumes was undertaken in the MATLAB-based CONN toolbox (version 19.c) (Whitfield-gabrieli and Nieto-castanon, 2012) to align images to an anatomical atlas and to reduce spatial and temporal artifacts. Structural scans were processed using the voxel-based quantification toolbox in SPM12 to generate MT and R1 maps (Weiskopf et al., 2013). The MT maps, which provide high contrast-to-noise ratio, were segmented to generate white matter, gray matter, and cerebrospinal fluid (CSF) maps using the “Segment” toolbox in SPM12 (Ashburner and Friston, 2005; Helms et al., 2009).

The first three volumes of each resting-state session were discarded. The remaining 253 functional volumes underwent realignment, unwarping, and slice-timing correction before undergoing rigid registration to the anatomical R1 scans and corresponding CSF and white-matter maps. Anatomical scans were then spatially normalized to Montreal Neurological Institute space (spatial resolution, 2 × 2 × 2 mm³); the resultant transformation was then applied to the functional data. Functional outlier detection was carried out using Artifact Detection Tools-based scrubbing¹. The functional volumes were then inspected and smoothing was applied to reduce potential spatial and temporal artifacts using an 8-mm full-width at half-maximum Gaussian kernel.

To reduce spurious sources of variance in the functional data, denoising was carried out using the component-based noise correction method (Behzadi et al., 2007). The temporal time series with estimated subject motion (average framewise displacement calculated from 3 rotation and 3 translation parameters) and the BOLD time series within the subject-specific white matter and CSF masks, were used as temporal covariates and removed from the BOLD data using linear regression. Ultra-low-frequency fluctuations in the resulting residual BOLD time series were removed using a high-pass filter at 0.0078 Hz.

An overview of the entire processing workflow is illustrated in Figure 1. In order to take a data-driven approach that would not be biased by the application of a specific structural atlas, independent component analysis (ICA) was used to extract spatially independent components (ICs) from the rs-fMRI data. The Minimum Description Length (MDL) algorithm was first used to determine optimal number of components to be extracted from the preprocessed rs-fMRI data of each patient; and the mean number of components across the cohort was calculated. Group principal component analysis (PCA) were performed first on a subject-specific level and then on a group level to reduce rs-fMRI data to this mean optimal number of components. Spatial group ICA was then performed by applying the Infomax algorithm to find the ICs; and stability was assessed by running the algorithm 20 times in the ICASSO toolbox (Himberg and Hyvärinen, 2003). Group ICA back reconstruction was then performed to generate aggregated spatial maps of the ICs for the cohort. All steps for group ICA were performed using the GIFT toolbox (version 3.0b) (Calhoun et al., 2001).

Manual classification of ICs as signal or noise remains the gold standard and as such, the resultant components were inspected independently by two of the authors (CW and CM) to identify and eliminate noise components according to previously described classification methods (Griffanti et al., 2017). Discrepancies were handled by a consensus decision between the two reviewers. Based on spatial location, the remaining components were then grouped into one of eight well-described large-scale functional networks: salience network (Sa); default mode network (DMN); frontoparietal network (FPN); sensorimotor network (SM); cerebellar network (Ce); visual network (Vi); language network (La); or auditory network (Au).

Using the extracted ICs as ROIs, overall static FC was calculated for each subject. Only the BOLD signal within each subject-specific gray-matter mask was included in these calculations. Bivariate correlation coefficients were calculated for the entire BOLD time series between each ROI pair (p < 0.05) to generate the FC matrices for each patient. Multivariate analysis of covariance (MANCOVA) was then used to assess and incorporate covariates into the group representation of overall static FC (Allen et al., 2011). Given the known FC changes that occur with age and gender, (Allen et al., 2011; Ferreira et al., 2015; Kim et al., 2017) age at the time of DBS implantation and gender served as a control covariates. Since images were acquired during the medication-ON state the preoperative UPDRS-III medication-ON score was also incorporated as a control covariate in order to account for disease severity. While some studies have also used LEDD as a covariate, we do not believe this to be a reasonable strategy because of the variable bioavailability of oral dopamine replacement therapy, (Di Stefano et al., 2009; Tuite, 2016) which yields different therapeutic effects across patients. Instead, the UPDRS-III medication-ON score best represents the clinical state resulting from optimal medical therapy and therefore, the effects of L-DOPA on FC. Lastly, the average framewise displacement calculated during image preprocessing was also used as a control covariate in order to minimize the potential contribution of patient movement at the time of rs-fMRI acquisition on FC (Johnstone et al., 2006). In summary, age, gender, preoperative UPDRS-III medication-ON score, and subject motion during scan acquisition were incorporated as covariates in the calculation of the overall static FC matrix for the entire cohort. These steps for FC analysis were performed using the MANCOVAN toolbox within the GIFT toolbox.

In preparation for graph theory analysis, the resulting group FC matrix underwent thresholding to remove spurious connections and create a sparse graph. In order to analyze an efficient network capable of specialized, distributed, and integrated information processing, the small-world topology of a graph should be maximized with this thresholding step (Bassett and Bullmore, 2006; Langer et al., 2013). Rather than selecting an absolute threshold value, the graphs were binarized by maintaining a percentage of the strongest connections. Such proportional thresholding has been suggested to more favorable because it avoids the generation of graphs with different network densities across the study cohort, which in turn may bias the resulting calculation of network statistics (van Wijk et al., 2010). At the same time, proportional thresholding remains susceptible to bias from differences in overall functional connectivity between subjects. Therefore, for each patient, the overall functional connectivity was calculated as the mean of all positive values of the connectivity matrix so that it could be incorporated as a control covariate into the partial correlation analysis (van den Heuvel et al., 2017).

While the optimal proportional threshold is typically associated with the graph with maximal small-worldness, (Bassett and Bullmore, 2006; Langer et al., 2013) to accommodate for differences in optimal thresholds between graphs from different patients, we calculated the area under the curve (AUC) of each graph theory metric of interest over a range of thresholds (Koshimori et al., 2016; Kim et al., 2017). The small-worldness across this threshold range was also calculated to ensure that the graph with maximal small-worldness was included in this analytic range. For each metric, a total of 25 undirected and unweighted graphs were calculated by applying thresholds ranging from 10–34% at increments of 1% (Achard and Bullmore, 2007). The AUC was then calculated for each graph theory metric.

Although such properties have been previously described in detail, (Bullmore and Sporns, 2009; Bullmore and Bassett, 2011) it is instructive to familiarize ourselves with the metrics employed in this study. Global network properties were assessed with global efficiency and assortativity. Global efficiency represents the average strength of connection between all pairs of nodes in the network and as such, represents the extent of overall network integration (Medaglia, 2017). Assortativity indicates the degree to which nodes with similar degrees of connectivity are connected to each other (e.g., highly connected nodes connected to other highly connected nodes and poorly connected nodes connected to other poorly connected nodes) (Newman, 2002). Networks that are more assortative are considered to be more resilient – as removal or malfunction of a node has a lesser effect on the function of the overall network (Medaglia, 2017). Overlapping metrics of integration (e.g., characteristic path length) and resilience (e.g., degree distribution) were not calculated to avoid redundancy.

Local network properties were assessed with clustering coefficients and betweenness centralities for each ROI. The clustering coefficient is a measure of the proportion of a node’s neighbors that are also directly connected to each other. This metric therefore serves as an indication of the local network segregation. Betweenness centrality highlights the relative importance of a node in the communication efficiency of local network (Zuo et al., 2012). It is calculated by determining the number of shortest paths of a network that involve that particular node. As such, nodes with high values of betweenness centrality participate in a large number of shortest paths and have an important role in the information transfer within a network (Nigro et al., 2016). Overlapping metrics of segregation (e.g., local efficiency) and centrality (e.g., closeness centrality) were not calculated to avoid redundancy.

These four complementary graph metrics were each calculated across the aforementioned threshold range to calculate the AUC for each metric. All graph analysis calculations were performed using the GRETNA toolbox (v2.0.0) in the MatLab environment (Wang et al., 2015).

To assess the relationship between imaging characteristics and STN-DBS response, Pearson partial correlations were performed. Since the age, motion, and preoperative UPDRS-III ON score were already incorporated as control covariates in the calculation of the FC matrix, these variables were not incorporated into the partial correlations. For statistical analysis of the AUC of the graph theory metrics, overall functional connectivity was also used as a control covariate for the reasons delineated above. In order to correct for multiple comparisons, p-values underwent FDR correction when evaluating correlations for FC and local graph theory metrics. Statistical analysis was performed using the RVAideMemoire (version 0.9-77) and stats (version 3.6.3) packages within R Statistical Software (version 3.6.3; R Foundation for Statistical Computing, Vienna, Austria).

Dynamic state classification was performed in a manner analogous to that described by the group from the University of Toronto using the GIFT toolbox (Kim et al., 2017; Díez-Cirarda et al., 2018). A sliding window of 22 TRs (59.4 s) with a 3 TR (8.1 s) Gaussian window and a step size of 1 TR (2.7 s) was applied to each subject’s rs-fMRI data, resulting in a total of 231 consecutive windows for dynamic analysis. For each window FC was calculated in a manner similar to that described for static FC; with the addition of L1 regularization to account for potential noise generated by the relatively short duration of the time windows. For each subject, 10 training repetitions were run in order to estimate the most appropriate regularization rate before application to the FC matrices.

Each of the 231 windows were then classified as a particular “brain state” based on the FC pattern by applying k-means clustering. To estimate the optimal number of brain states for classification, silhouette analysis was performed while varying the number of clusters between 2 and 10. The algorithm was run 100 times to reduce bias from random selection of cluster centroids. Bayesian information criterion (BIC) based on the Euclidean distance between clusters was then used to validate the result of the k-means clustering. In this manner, all 231 FC matrices for each subject could be appropriately classified into a particular brain state. The median FC values for all matrices classified as a particular brain state were calculated to generate a single d-FC matrix for that brain state. To better understand the significance of these states, the variance of FC across all subjects was also calculated; and compared to the variance in static FC. In addition, temporal properties of dwell times and number of state transitions were also calculated for further statistical analysis.

For each d-FC matrix, graph metrics were calculated in the manner described above for static analysis. Once again, the AUC for the four graph theory metrics of interest were calculated with thresholds from 10–34%. It is again important to note that because small-worldness may be optimized at different thresholds for graphs across different subjects, this approach also allows for a uniform comparison of metrics across different graphs representing different brain states.

Partial correlations were performed to assess the relationship between the temporal properties of d-FC (e.g., dwell times, number of state transitions) and clinical responses to STN-DBS. Analysis of d-FC itself and of the calculated graph metrics were performed for each state in the same manner described above for static analysis.

Application of the MDL algorithm to the preprocessed rs-fMRI data yielded a mean optimal number of 66.579 ± 11.801 independent components. As such, data reduction with PCA reduced data first to 101 principal components based on subject-specific data, before arriving at 67 group independent components. After ICA and back reconstruction, the total number of ICs was reduced to 25 following manual classification. The 25 ICs were classified into 8 networks (Figure 2 and Table 2).

The overall static FC for the entire cohort is represented as the leftmost matrix in Figure 3. ROIs generally appear strongly interconnected – particularly ROIs within the DMN, SM, FPN, and Ce networks. There were no statistically significant correlations between FC and response to STN-DBS at one year after DBS implantation. Further analysis with graph theory metrics of global efficiency, assortativity, clustering coefficients, and betweenness centralities also failed to demonstrate any significant relationships with STN-DBS.

Using both silhouette analysis and BIC for k-means clustering, categorization of d-FC into one of two brain states was determined to be optimal (Supplementary Figure 1).

Subjects generally spent more time overall in State 2 (67.03 ± 32.64%) than in State 1 (32.97 ± 32.64%); and also spent longer contiguous blocks of time in State 2 (4.62 ± 4.07 min) than in State 1 (1.00 ± 1.02 min). Six of the 18 patients spent more time in State 1 than in State 2. Across the cohort, patients transitioned between states an average of 4 ± 2.77 times. Time spent in a particular state, dwell time, and number of transitions were not significantly correlated with response to STN-DBS.

The representative FC matrices for the entire cohort for each brain state are illustrated on the right-hand panels of Figure 3. While the d-FC State 1 appears similar to the overall static FC and relatively hyperconnected, State 2 represents a relatively hypoconnected state, with stronger connections only remaining in portions of the SM, DMN, and Vi networks.

Of note, variance in FC across subjects decreased by an order of magnitude for both d-FC states (Figure 4). While mean variance for the total timecourse was 0.204 (range 0 – 0.541), mean variance for State 1 was 0.022 (range 0 – 0.074) and mean variance for State 2 was 0.016 (range 0 – 0.053).

There were no statistically significant correlations between d-FC measures and response to STN-DBS at one year after DBS implantation. Further analysis with graph theory metrics demonstrated significant relationships with STN-DBS only in the hypoconnected state (State 2). Specifically, response to STN-DBS was negatively correlated with network assortativity (r = −0.6907 [95%CI = −0.8937 −0.2528], p = 0.006) (Figure 5). No other graph metrics were significantly correlated with STN-DBS.

For patient comfort and to reduce motion artifact, all patients were on optimal L-DOPA medication at the time of the preoperative MRI scan. As such, we were not able to directly assess the effects of L-DOPA and DBS on BOLD signal. At the same time, since our aim was to evaluate the potential of rs-fMRI to predict treatment to DBS, we only used data that could be reasonably obtained only prior to surgery.

While the relatively small sample size may be considered a limitation of this work, the 18 subjects enrolled in the current study is comparable to other rs-fMRI studies involving PD patients (Tahmasian et al., 2015). Furthermore, despite this potential limitation, the identified correlations are strong and certainly warrant further investigation. Another limitation of this work is the relative heterogeneity of the study population. Since there exists no standard diagnostic criteria for characterizing advanced PD, (Fasano et al., 2019) there still exists variability in terms of age, disease severity, and motor outcomes. Nevertheless, we did control for such factors in an attempt to perform the most robust analysis possible with the existing data. Ultimately, we believe that the heterogeneity of our cohort benefits this study, as the results are more generalizable to the population of interest – namely those that may be potential candidates for DBS.

It is also important to note that no ICs were identified within the basal ganglia. While the reason for this deficiency remains unclear, the data-driven approach we have taken in our analysis prevents us from forcing the inclusion of this region. Even upon careful review of the 67 original ICs, there was no IC with meaningful signal within the basal ganglia. We must therefore keep this in mind in the interpretation of our results, as the importance of the basal ganglia cannot be refuted in the study of PD. In the current study, we are unable to make any conclusions about their role in predicting STN-DBS response.

The lack of basal ganglia ICs is one reason we are unable to make direct comparisons with our prior work on L-DOPA response in an overlapping cohort of patients (Akram et al., 2017). Given the overlapping cohort used in the present study and our prior work on FC in predicting L-DOPA response, it would be useful to draw comparisons between these two analyses – particularly given the known relationship between L-DOPA and DBS response. Significant differences in methodology, however, also prevent us from making such comparisons. Our prior work, employed a hypothesis-driven approach in which a narrow scope of analysis was performed solely on predefined regions of interest defined on the FSL Harvard-Oxford cortical structural atlas, the ATAG subthalamic nucleus atlas, and the ATAG MNI04 BG atlas. Static functional connectivity was then analyzed to investigate relationships between FC and L-DOPA response. The current study, however, takes a data-driven approach using ICA-defined ROIs, which generates different regions for analysis than with our prior work. Furthermore, instead of seed-to-seed analysis in which individual connections are being analyzed, we have taken a whole-brain network approach with graph theory analysis. We have also employed dynamic resting-state analysis to further the sensitivity of the analysis by incorporating the investigation of different brain states. With that said, if independent components had been identified within the basal ganglia, we could potentially have made direct comparisons with this prior work.

Given the strict criteria implemented for statistical significance, it is certainly possible that there are other relevant relationships that have been overlooked due to the sample size of our cohort and type II error. Future endeavors to replicate this study or with larger cohorts capable of providing more statistical power will provide a better understanding of these relationships.