Alejandro N. Santos,Ferath Kherif,Lester Melie-Garcia,Antoine Lutti,Alessio Chiappini,Laurèl Rauschenbach,Thiemo F. Dinger,Christoph Riess,Amir El Rahal,Marvin Darkwah Oppong,Ulrich Sure,Philipp Dammann,Bogdan Draganski
Parkinson’s disease may disrupt overlapping subthalamic nucleus and pallidal motor networks
Parkinson’s disease,Multi-parameter mapping,Structural covariance,Magnetic resonance imaging,Voxel-based morphometry,Voxel-based quantification
Highlights • We use MRI to obtain whole-brain maps sensitive to myelin, iron and tissue water. • We show overlapping subthalamic nucleus and internal pallidum covariance patterns. • The patterns of tissue property covariance are confined to the motor system. • We show differential covariance patterns in patients with Parkinson’s disease.
Parkinson’s disease may disrupt overlapping subthalamic nucleus and pallidal motor networks

• We use MRI to obtain whole-brain maps sensitive to myelin, iron and tissue water. • We show overlapping subthalamic nucleus and internal pallidum covariance patterns. • The patterns of tissue property covariance are confined to the motor system. • We show differential covariance patterns in patients with Parkinson’s disease.

Given the critical role of the basal ganglia and thalamus in modulating motor and non-motor cortico-subcortical networks in movement disorders, there is a pressing need to delineate the connectivity patterns determining their unique contributions to clinical outcomes. The subthalamic nucleus (STN) and globus pallidus pars interna (GPi) are established targets for deep brain stimulation (DBS) in patients with idiopathic Parkinson’s disease (PD) (Kleiner-Fisman et al., 2006, Volkmann et al., 1998). Despite significant progress in the field, we are still falling short of understanding the differential effects of DBS linked to a particular target structure, which hampers the accurate prediction of clinical outcome and the prevention of stimulation-associated adverse effects (Lachenmayer et al., 2021). Most recent meta-analytical evidence (Lachenmayer et al., 2021), covering 14 years after the first STN vs GPi meta-analysis and including over 2000 DBS-treated patients with PD (Kleiner-Fisman et al., 2006), supports the notion of overall non-superiority between targets on the background of significant disparity of the available pre- and post-surgery data. With respect to therapy efficacy, the results point towards better motor outcome after STN DBS, whilst the fair comparison for other treatment end-points - non-motor symptoms and quality-of-life, was not possible due to the heterogeneity and scarcity of data, particularly for GPi DBS. The anatomical embedding of the STN and GPi in the cortico-subcortical circuitry favours their partially overlapping functional profiles (Shen et al., 2020, Younce et al., 2021, Zhang et al., 2021, Haynes and Haber, 2013). Although distinct structures, STN and GPi are strongly interconnected via the bidirectional subthalamic fasciculus. Through the hyper-direct pathway, the dorso-lateral STN receives premotor and motor cortical projections (Nambu et al., 1996), whilst the GPi projects via the ventro-lateral thalamus to the same premotor and cortical motor areas (Nakano et al., 2000). With the advent of non-invasive magnetic resonance imaging (MRI), diffusion-weighted tractography has taken the lead in investigating the anatomical connectivity of brain structures over the traditional tract-tracing techniques in the animal model. Currently, diffusion-based tractography is increasingly used for enhancing the precision of individuals’ DBS electrode placement (Liebrand et al., 2019). However, the downside of the available tractography methods is the presumption of representing only monosynaptic connections and their reliance on relatively low-resolution data – i.e., 1.5–2 mm resolution on a clinical 3Tesla scanner. These, together with the inability to differentiate between dominant and small bundles’ contributions result in a high rate of false positives (Maier-Hein et al., 2017) that renders the technique unsuitable for in vivo discovery of new pathways. The investigation of structural covariance overcomes part of the diffusion-based tractography shortcomings by allowing for whole-brain investigation at the voxel/vertex or region-of-interest level at the usual 1 mm spatial resolution (Evans, 2013). Independent of the applied analytical method and features extracted from the MRI data at hand, the neurobiological interpretation of the obtained covariance networks also includes polysynaptic connections that go beyond the delineation of diffusion-based direct fiber pathways (Seeley et al., 2009). Initially established for studying covariance networks of cortical thickness or regional grey matter volume, this approach was extended by our group for the analysis of MRI measures indicative for microscopic brain tissue properties – i.e., myelin and iron, across cortical and subcortical structures (Accolla et al., 2014, Melie‐Garcia et al., 2018, Weitnauer et al., 2021). We sought to tackle the question about STN and GPi connectivity raised recently in a functional connectivity study (Sobesky et al., 2022), which tends to generalize across basal ganglia disorders (Haber et al., 2020, Haber et al., 2021). The supposition here is that the cortical and subcortical areas constituting a functional network will also share a similar anatomical fingerprint. We extend the existing morphometry description of the regional volume and cortical thickness to MR contrast estimates sensitive to myelin, iron and tissue water content (Accolla et al., 2014, Melie‐Garcia et al., 2018, Weitnauer et al., 2021). First, we compare the STN and GPi structural covariance patterns in a normative large-scale community-dwelling population (n = 1184). Next, we look for structural covariance differences in a cohort of patients with PD (n = 32) not only to provide face validity of our analytical framework but also to look for informative patterns associated with the view on PD as a disconnection disorder. Finally, we probe the robustness of the obtained patterns in a reduced size healthy controls’ cohort (n = 32).

We use data from the community-dwelling population (n = 1184, 629 females, mean age: 59,6 years, standard deviation ± 7,29 years) collected within the CoLaus|PsyColaus longitudinal study with more than 15 years of observation every five years. In addition, we study patients with Parkinson’s disease (PD) recruited at the local Movement Disorders speciality clinic (n = 32, 13 females, mean age: 64,4 years, standard deviation ± 6,41 years) that was partially presented in a previous publication (Jastrzębowska et al., 2019). For validation, we created a reduced size sex- and age-matched random subsample of the community-dwelling cohort (n = 32, 15 females, mean age: 65,9 years, standard deviation ± 3,92 years). For details on demographic and clinical characteristics, see Table 1. The local ethics committee approved the use of the anonymized data for further analysis beyond the initially intended study. Study data are available upon reasonable request.

The MRI data were acquired at a 3 T whole-body MRI system (Prisma, Siemens Healthcare, Erlangen, Germany) using a 64-channel RF head receive coil, and RF body transmit coil. The whole-brain relaxometry protocol comprised 3D multi-echo FLASH datasets with predominantly proton density weighting (PDw; repetition time TR = 23.7 ms, flip angle α = 6°), T1 weighting (T1w; TR/α = 18.7 ms/20°), and magnetization transfer weighting (MTw; TR/α = 23.7 ms/6°) contrast according to the previously published protocol (Melie‐Garcia et al., 2018). We acquired multiple gradient echoes with alternating readout polarity at six equidistant echo times (TE) between 2.34 and 14.04 ms for the MTw acquisitions and at 8 equidistant TE between 2.34 ms and 18.72 ms for the T1w and PDw acquisition. The Image resolution was 1x1x1mm voxel size, FOV 256 × 240 × 176 mm, matrix 256 × 240 × 176, GRAPPA factor 2 in phase-encoding (PE) direction, 6/8 partial Fourier in partition direction, non-selective RF excitation. For the correction of inhomogeneous radiofrequency (RF) excitation effects on the qMRI estimates, B1 + -mapping data was acquired using a 3D echoplanar imaging spin-echo (SE)/stimulated echo (STE) method with a FOV 256x192x192mm, matrix 64x48x48, TR = 500 ms (Lutti et al., 2010, Lutti et al., 2012). We corrected the B1 + maps for EPI image distortions and off-resonance effects using a standard B0 map according to the published protocol (Lutti et al., 2010, Lutti et al., 2012). All MRI data for the CoLaus|PsyCoLaus and the PD cohorts were acquired at the very same MRI scanner and under identical Syngo software VD13 version according to the same imaging protocol.

Given the effects of head motion on image quality and analysis particularly in the context of large-scale studies and in patients with movement disorders, i.e. Parkinson’s disease, we quantitatively estimate MRI data quality in our study. We assess the quality of the MRI data using the validated Motion Degradation Index (MDI) (Castella et al., 2018). Consistently with (Lutti et al., 2022), MDI values from the MT-weighted images are higher than those from the PD- and T1-weighted images by ∼ 1 s-1* due to the lower number of echo images and lower signal-to-noise of the raw data. The distribution of the MDI values across both cohorts is provided as Supplementary material (supplementary Figure 2). Datasets with MDI values below 6 s−1 for the PD- and T1-weighted images, and 7 s−1 for the MT-weighted images are empirically considered of sufficient quality. Extension of this method to such covariance analyses is currently under investigation. For completion, we also provide values of the rigid body transformation parameters for co-registration of the MT-weighted and T1-weighted images to the reference PD-weighted images, that represent metrics of head motion between image volumes (Tabelow et al., 2019) (supplementary Figure 3).

Image processing was performed with the freely available Statistical Parametric Mapping software (SPM12; Wellcome Centre for Human Neuroimaging, London, UK, https://www.fil.ion.ucl.ac.uk/spm/software) and inhouse routines, running under Matlab2019a (Mathworks, Sherborn, MA, USA). Multi-parameter maps calculation: Maps of magnetization transfer (MT) saturation and effective proton density (PD*) as described in (Helms et al., 2008, Helms et al., 2008), from the raw MR images averaged across echo times (Tabelow et al., 2019, Draganski et al., 2011). Maps of the effective transverse relaxation rate (R2* = 1/T2*) were computed from the raw MR images pooled across all contrasts using the ordinary-least-squares approach described in (Weiskopf et al., 2014). Following this, we automatically classified the brain images into grey matter, white matter and cerebro-spinal fluid (CSF) using the MT and PD* within the multi-channel”unified segmentation” framework of SPM12 with enhanced tissue priors for optimal delineation of basal ganglia (Lorio et al., 2016). The individual grey and white matter tissue classes were used for diffeomorphic registration to standard Montreal Neurological Institute (MNI) space (Ashburner, 2007). We applied the spatial registration parameters to the parameter maps following the weighted averaging procedure of the voxel-based quantification (VBQ) framework (Draganski et al., 2011). Finally, we apply an isotropic Gaussian kernel of 8 mm full-width-at-half-maximum for spatial smoothing before statistical analysis. Motor STN and GPi seeds: Aiming to test the hypothesis of structural covariance within motor networks, we used a previously published definition of the motor part of the STN (mSTN) and GPi based on probabilistic diffusion-based tractography (Accolla et al., 2014). We then extracted from the volume and multi-parameter maps registered to the standard MNI space the mSTN and GPi local averages of PD*, MT saturation, R1, R2* values indicative for tissue water, myelin and iron content that were used for whole-brain covariance analysis.

For the whole-brain structural covariance analyses, we used a linear regression model as implemented in the General Linear Model framework of SPM. Age, gender, and total intracranial volume (TIV) were included as additional variables in the design matrix to control for global effects on the brain structure. Voxel-based two-tailed T-statistics were calculated to detect GPi and mSTN covariance patterns for each parameter map and each seed structure separately. Statistical thresholds were applied at p < 0.05 after family-wise error (FWE) correction for multiple comparisons over the entire GM or WM search volume. Trends were assessed using an auxiliary uncorrected voxel threshold of p < 0.001. The use of parametric tests for our analyses is motivated by the large sample theorem, which states that as the sample size increases, the sample mean approaches the population mean and the sample distribution becomes more normal, regardless of the population distribution. The law applies to both parametric and non-parametric tests if the sample size is large enough. Non-parametric approaches make no assumptions about the underlying distribution, but they also require a larger sample size. Non-parametric tests can be less efficient and powerful than parametric tests if the sample size is large. In the case of large samples, parametric methods are more efficient because they can build on the assumed structure to obtain more accurate and robust estimates (Flandin and Friston, 2019). In the case of functional MRI (fMRI) data, the complexity of the noise and the assumptions about the shape of the variance matrix of the noise can be problematic for parametric tests (and non-parametric tests), which can result in systematic bias (Winkler et al., 2014). In brain anatomy studies, the structure of the noise is assumed to be independent across individuals, so we benefit from the law of large numbers (Salmond et al., 2002).

Volume estimates: We observed largely overlapping structural covariance patterns of the mSTN and GPi, including the basal ganglia, the primary motor cortex and insula bilaterally (all results at pFWE < 0.05; see Fig. 1 and supplementary Table 1-2). The covariance between the local mSTN and GPi volume and white matter comprised all three components of the superior longitudinal fasciculus (all results at pFWE < 0.05; see Fig. 1 and supplementary Table 1-2). In the reduced size cohort, we could replicate only the subcortical and the midline cortical findings (supplementary Figure 3). Magnetization transfer saturation MTsat: The MT saturation in the mSTN and GPi covaried with basal ganglia, thalamus, hippocampus, entorhinal cortex, insula and primary motor cortex bilaterally (all results at pFWE < 0.05; see Fig. 2 and supplementary Table 1-2). Across the white matter, MT saturation in the mSTN and GPi showed significant covariance with the temporal pole, the posterior limb of the internal capsule and the white matter adjacent to the primary sensorimotor cortex (all results at pFWE < 0.05; see Fig. 2 and supplementary Table 1-2). Similarly, in the reduced size cohort, we could replicate only the subcortical and the midline cortical findings (supplementary Figure 4). Effective proton density PD*: Local PD* averages in the GPi covaried with the PD* content in putamen, external pallidum, insula and cingulate gyrus in both hemispheres, whereas the mSTN covaried with thalamus, red nucleus and substantia nigra (all results at pFWE < 0.05; see Fig. 2 and supplementary Table 1-2). Transverse relaxation rate R2*: The R2* averages in the mSTN and GPi showed a more spatially restricted pattern confined to the basal ganglia, the cerebellar dentate nucleus and the surrounding white matter R2* content (all results at pFWE < 0.05; see Fig. 3 and supplementary Table 1-2).

Volume estimates: Like the findings in healthy controls; we observed an overlapping structural covariance pattern of the mSTN and GPi that included the basal ganglia and surrounding white matter without further extending to the cortex (all results at pFWE < 0.05; see Fig. 1 and supplementary Table 3-4). Magnetization transfer saturation MTsat: We report a covariance between mSTNs MT saturation with the MT saturation in basal ganglia, thalamus, hippocampus, entorhinal cortex, insula and planum temporale bilaterally (all results at pFWE < 0.05; see Fig. 2 and supplementary Table 3-4). The GPi did not show any significant covariance results. Similarly, we did not observe any covariance between the MT saturation of mSTN and GPi across the white matter. Effective proton density PD*: Conversely to the PD* covariance pattern confined to subcortical areas in healthy controls, PD patients show an extended anatomical pattern including not only basal ganglia and thalamus, but also the cingulate gyrus, insula, precuneus, cerebellar and mesial temporal lobe structures for both mSTN and GPi (all results at pFWE < 0.05; see Fig. 2 and supplementary Table 3-4). Transverse relaxation rate R2*: R2* averages in the mSTN and GPi covaried with the basal ganglia, the cerebellar dentate nucleus (GPi only) and the surrounding white matter R2* content (all results at pFWE < 0.05; see Fig. 3 and supplementary Table 3-4).

Our structural covariance study leveraged the statistical power of large-scale data (n greater than 1000) to provide empirical evidence about the overlapping anatomical connectivity patterns of both DBS targets for treatment of PD – the GPi and the mSTN. We test the subcortical covariance patterns in a separate clinical cohort of PD patients and demonstrate a differential anatomical profile of cortical regions. While GPi’s and mSTNs volume and myelin content in healthy controls covary with the corresponding tissue properties of the primary motor cortex, PD patients seem to lose it and show covariance of increased unbound water content in these very same GPi and mSTN related motor networks. The partial validation of the cortical covariance patterns in a reduced size cohort of healthy controls limited to midline motor areas only, calls for a cautious interpretation of the obtained findings. Our main finding is the discovery of almost identical structural covariance patterns of both DBS targets for the treatment of PDs motor fluctuations – the GPi and mSTN across the brain. Similar to our previous findings, the delineated networks comprise subcortical areas representing key nodes of the motor circuitry (Accolla et al., 2014). The novelty is the demonstration of covariance with cortical primary motor areas – for both regional volume and indices of myelin content in the GPi and mSTN. These results align with relaxometry and myelin water fraction differences in patients with PD (Dean et al., 2016) and more recent work on the diffusion imaging-derived PD connectome weighted by local R1 values (Boshkovski et al., 2022). The similar covariance pattern of GPi and mSTN confirms the notion of shared tissue properties across cortico-subcortical networks that go beyond monosynaptic connections – e.g. the hyperdirect pathway between the cortex and mSTN. Finally, our findings corroborate the presumed overlap of functional connectivity patterns between the GPi and mSTN as recently demonstrated on the empirical basis of predictability of STN-DBS outcome by GPi connectivity and vice versa (Sobesky et al., 2022). Conversely, the R2* covariance maps indicative of tissue iron content remained confined to subcortical and cerebellar structures for both DBS targets. We attribute this finding to the specific distribution of dopaminergic receptors, intra- and extra-striatal dopamine across the basal ganglia (Smith and Kieval, 2000, Smith and Villalba, 2008). The presumption here is that the involvement of iron in the enzymatic machinery for dopamine synthesis and reuptake along nigro-striatal, nigro-pallidal and nigro-subthalamic projections, additionally to its storage in nigral neuromelanin and ferritin renders the basal ganglia as susceptible to iron-related MR contrast behaviour compared with the cerebral cortex (Ward et al., 2014). The second finding of our study is the loss of structural covariance of GPi’s and mSTNs volume and MT saturation with primary motor areas in PD patients. This is contrasted by the emergence of widespread PD* covariance with cortical and subcortical structures that is much less spatially extended in healthy controls. Under the assumption that neurodegeneration is related to increase in MR visible water due to cell loss, thus represented by increase in PD*, we interpret the concomitant loss of covariance with primary motor cortex and increased covariance as disruption of motor networks in PD. The differential effects on volume, MT saturation as index of myelin content and PD* lend us some confidence in the validity of our results. We interpret our results of loss of volume and MT saturation covariance between the GPi, mSTN and the superior longitudinal fasciculus in the context of neurodegenerative processes governing PD progression. Similar to the notion of shared mechanisms underlying the phenomenon of structural covariance in grey matter (Evans, 2013), atrophy or loss of myelin in the GPi or the mSTN are associated with volume decrease in long-range fronto-parietal connections, as shown previously (Dadar et al., 2020). The covariance across brain tissue compartments also provides neurobiological face validity to the proposed analysis within networks rather than specific tissue types. Despite the sufficient statistical power, we would like to acknowledge certain limitations of our study. First, one can argue that the absence of volume and MT saturation structural covariance between the GPi, mSTN and primary motor cortex in PD is merely the issue of statistical power rather than true loss due to neurodegeneration. This assumption is supported by the fact that in the reduced size cohort we could validate the covariance patterns only in midline motor areas. The differential impact of PD on atrophy/myelin vs. effective proton density covariance could, however, point towards true biological effects. Additionally, even lowering the statistical threshold to uncorrected levels did not show any trends for covariance with primary motor areas in PD patients. Second, we also acknowledge potential bias in the PD* and R2* maps in iron-rich structures prone to susceptibility artefacts due to violation of the assumption of the mono-exponential decay model (Weiskopf et al., 2013). Given the fact that we do not expect a differential impact of this potential bias in healthy controls and PD patients, we estimate its effect as negligible. Finally, given the still ongoing work on implementing statistical methods for structural covariance analysis of multiparameter maps, we acknowledge the lack of correction for the effects of head motion on the statistical validity as recently suggested (Lutti et al., 2022). In summary, we investigate the structural covariance patterns of the two most widely used DBS targets in Parkinson’s disease – the GPi and mSTN, across multiple available brain anatomy characteristics ranging from morphometry via parameters indicative for tissue myelin, iron, and unbound water content. We demonstrate the spatial pattern of covariance in a large-scale community dwelling cohort that spans over cortical and subcortical nodes of motor networks. Confirmatory to recent studies about largely overlapping functional connectivity between the GPi and mSTN, we find a similar structural covariance pattern. The analysis in a clinical cohort of patients with PD demonstrates a loss of the myelin and volume-defined covariance with primary motor cortex paralleled by emerging covariance patterns of effective proton density related to the process of neurodegeneration. However, the fact that we replicate only partially the cortical covariance patterns in a reduced size cohort of healthy controls, motivates the cautious interpretation of the obtained findings as a signature of PD-associated disruption of motor circuits.

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.