Changes in functional organization and white matter integrity in the connectome in Parkinson's disease.

Parkinson's disease (PD) leads to dysfunction in multiple cortico-striatal circuits. The neurodegeneration has also been associated with impaired white matter integrity. This structural and functional "disconnection" in PD needs further characterization. We investigated the structural and functional organization of the PD whole brain connectome consisting of 200 nodes using diffusion tensor imaging and resting-state functional MRI, respectively. Data from 20 non-demented PD patients on dopaminergic medication and 20 matched controls were analyzed using graph theory-based methods. We focused on node strength, clustering coefficient, and local efficiency as measures of local network properties; and network modularity as a measure of information flow. PD patients showed reduced white matter connectivity in frontoparietal-striatal nodes compared to controls, but no change in modular organization of the white matter tracts. PD group also showed reduction in functional local network metrics in many nodes distributed across the connectome. There was also decreased functional modularity in the core cognitive networks including the default mode and dorsal attention networks, and sensorimotor network, as well as a lack of modular distinction in the orbitofrontal and basal ganglia nodes in the PD group compared to controls. Our results suggest that despite subtle white matter connectivity changes, the overall structural organization of the PD connectome remains robust at relatively early disease stages. However, there is a breakdown in the functional modular organization of the PD connectome.

Introduction

Parkinson's disease (PD) is a neurodegenerative disorder characterized by dopaminergic neuronal loss in the substantia nigra pars compacta (Braak et al., 2003). Dopamine deficiency leads to dysfunction in multiple, topographically organized cortico-striatal circuits, consequently giving rise to motor symptoms as well as cognitive and emotional problems. The pathology in PD is not limited to the neuronal cell bodies. Studies using animal and cell culture models of PD have demonstrated that neurodegeneration is also associated with axonopathy and synaptic dysfunction, and impairs white matter integrity (Burke and O'Malley, 2013; O'Malley, 2010; Tagliaferro et al., 2015).

PD-related deficits are complex in nature. The neuroanatomical correlates of these wide-ranging deficits have been characterized in task-based neuroimaging studies (for review Hanganu et al., 2015; Schendan et al., 2013; Tinaz et al., 2008). Several studies have also investigated the alterations in functional connectivity patterns in PD using resting-state functional MRI (rs-fMRI). Many of these studies focused on abnormal functional connectivity in cortico-striatal circuits (Hacker et al., 2012; Helmich et al., 2010) or in specific networks (e.g., default mode, dorsal attention, salience) using various methods such as seed-based functional connectivity or independent component analysis (Baggio et al., 2015; Putcha et al., 2015; Tessitore et al., 2012). Neuroimaging studies using diffusion tensor imaging (DTI) have shown disrupted white matter integrity in major tracts in various stages of PD. Furthermore, this disruption was found to correlate with cognitive impairment (Hattori et al., 2012; Melzer et al., 2013).

More recently, complex network theories and tools have contributed substantially to our understanding of the structure and function of the brain as a large-scale neural network, i.e., connectome. Within this framework, pathology is approached as a phenomenon emerging from the abnormal connections and interactions between distributed brain regions rather than being the result of focal lesions. Network research using neuroimaging data has led to important insights into the brain pathology in several neuropsychiatric conditions including Alzheimer's disease, schizophrenia, traumatic brain injury, epilepsy, and multiple sclerosis (Filippi et al., 2013; Stam, 2014). Relatively few studies have employed large-scale graph theory-based network methods using rs-fMRI in PD. These studies reported abnormalities in global and local network properties within the studied networks which also correlated with disease severity and cognitive deficits (Berman et al., 2016; Göttlich et al., 2013; Lebedev et al., 2014; Skidmore et al., 2011; Tinaz et al., 2016).

Taken together, evidence from neuropathological and neuroimaging studies suggests that PD can be conceptualized as a structural and functional “disconnection” syndrome. Yet, the nature of this disconnection warrants better characterization. Results of the rs-fMRI network studies in PD have been variable due to methodological factors and cohort characteristics, and the whole connectome was not examined in all studies. In addition, to our knowledge, the white matter integrity of the whole connectome in PD has not been investigated using graph theory-based network methods.

Our main goal in this study was to investigate the white matter integrity and functional reorganization of the connectome in non-demented, medicated PD patients using DTI and rs-fMRI, respectively. We also examined whether white matter and functional connectivity changes coincided spatially within the connectome. To this end, we employed graph theory-based methods specifically focusing on local network metrics including node strength, clustering coefficient, and local efficiency; as well as network modularity.

We hypothesized that there would be nodal connectivity changes and breakdown in the modular organization of the structural and functional networks in PD patients compared to matched controls.

Materials and methods

Subjects

The following exclusion criteria applied to all participants: The presence of any neurological or psychiatric disorder (other than PD), abnormality in neurological examinations or in routine clinical MRI scans, and active alcohol or illicit drug use.

Patients were referred from the National Institutes of Health Parkinson's Disease Clinic. Healthy volunteers (HVs) were recruited from the local community. Twenty-four patients with idiopathic PD (age range 44–75, 12 females) and 24 matched healthy volunteers (HV) (age range 43–71, 11 females) participated in the study after giving written informed consent in accordance with the Combined Neuroscience Institutional Review Board of the National Institutes of Health.

The PD diagnosis was established according to the UK Parkinson's Disease Society Brain Bank Clinical Diagnosis Criteria (Hughes et al., 2001). All patients had bradykinesia and at least one of the following impairments: rigidity, resting tremor, or postural instability. Patients were assessed using the Unified Parkinson's Disease Rating Scale (UPDRS) (Fahn and Elton, 1987) and the Hoehn and Yahr (H&Y) scale (Hoehn and Yahr, 1967). The Mini Mental State Examination (MMSE) (cut-off score < 26) (Dubois et al., 2007) and Montreal Cognitive Assessment test (MoCA) (Nasreddine et al., 2005) (cut-off score < 19) (Hoops et al., 2009) were also administered to all participants to screen for dementia.

Scanning procedure

All patients were clinically evaluated using the UPDRS part III (motor exam) first off any dopaminergic medication > 12 h in the morning. Patients took their regular dose of dopaminergic medications immediately before they were placed in the MRI scanner. Structural scans were collected first, all of which lasted about an hour. Subsequently, in the same session, the rs-fMRI scans were collected within 60–70 min following medication intake. The timing of the rs-fMRI scans was planned to allow to reach the clinically defined “on” state which is usually observed in one hour following medication intake. The clinical assessment during the “on” state was performed immediately upon completion of scanning, within 80 min of the dopaminergic medication intake. The same scanning order was used for HVs.

Final cohort

Based on our criteria, three patients had excessive head motion (average motion ≥ 0.4 mm) in MRI scans and were excluded. In addition, one patient was excluded because of technical problems in MRI data collection and one HV was excluded due to a previously unkown lacunar stroke observed in the clinical MRI sequences. Three matched HVs were also excluded to match the PD group. In the end, a total of 20 patients (age range 44–75, average age 62.5 ± 6.9, 11 females) and 20 HVs (age range 43–71, average age 61.9 ± 6.6, 11 females) were included in the imaging analyses.

Image acquisition

All images were collected in a GE 750 3 T scanner using a 32-channel head coil in the Nuclear Magnetic Resonance Center at the NIH.

The following imaging parameters were used:

T1-weighted anatomical images: Magnetization prepared rapid acquisition gradient echo (MPRAGE), 3D inversion recovery, TR: 7.664 ms, TE: 3.42 ms, TI: 425 ms, slice thickness: 1 mm, 1 × 1 mm in-plane resolution, percent phase FoV: 100, flip angle: 7, matrix size: 256 × 256.

T2-weighted diffusion images: Fat-saturated, 62 images, TR: 7.500 ms, TE: 100.74 ms, slice thickness: 2.5 mm, 0.9375 mm × 0.9375 mm in-plane resolution, FoV: 240, percent phase FoV: 80, flip angle: 90, matrix size 256 × 256, maximum B value 1100, 70 diffusion directions and 10 non-diffusion weighted images.

Functional MRI data: Echoplanar images (EPI), TR: 2000 ms, TE: 27 ms, 40 slices with thickness: 3 mm, 2.5 × 2.5 mm in-plane resolution, axial orientation, FoV: 240, percent phase FoV: 100, flip angle: 70, matrix size: 96 × 96. FMRI data were collected during rest with eyes closed for 5 min.

Image processing

Pre-processing

DTI data

We used the same preprocessing and probabilistic tractography methods as detailed in Lauro et al. (2016). The T2-weighted images were mid-sagittally aligned and the anterior and posterior commissure (AC and PC) landmarks were manually defined on the T2-weighted volume using the Medical Image Processing, Analysis and Visualization (MIPAV) software package (McAuliffe et al., 2001). A rigid body Talairach transformation was then applied resulting in AC–PC aligned volumes. The T2-weighted images served as the co-registration target for the DWI volumes and the T1-weighted volume.

DWI volumes were used to estimate diffusion tensors and maps of their associated DTI parameters including the principal directions of diffusion (i.e. eigenvectors), eigenvalues, fractional anisotropy (FA), and mean diffusivity which were used to quantitatively compare brain tissue properties. FA and principal diffusion directions were used for estimating white matter fiber trajectories which estimate the proxy location of the white matter and the orientation of propagating tracts, respectively, and their confidence intervals were calculated for the implementation of probabilistic tractography.

Preprocessing of the DWI volumes was performed in TORTOISE (Pierpaoli et al., 2010). The DWI volumes were first motion-, eddy-, and EPI distortion-corrected, co-registered to the AC–PC aligned T2-weighted volume, and resampled to 1.5 mm isotropic voxels using standard settings in the ‘DIFF_ PREP’ tool.

The FA and directionally encoded color data were subsequently computed using ‘DIFF_CALC’, and resultant datasets were visually quality checked and exported to nifti format.

The gradient direction data obtained from DWI DICOM images were translated into AFNI-readable format using the 1dDW_Grad_o_Mat command with scanner-specific formatting. The preprocessed DWI volumes were used to estimate the diffusion tensors using AFNI's ‘3dDWItoDT’ command with nonlinear fitting.

We used the Functional and Tractographic Connectivity Analysis Toolbox (FATCAT) (Taylor and Saad, 2013) for tractography analyses. In FATCAT, both deterministic and probabilistic tractography use the algorithm called “fiber assessment by continuous tracking including diagonals (FACTID)”. This algorithm has been validated in a series of tests on both human data and a standard phantom (Taylor et al., 2012). Whole brain deterministic tractography was implemented using the ‘3dTrackID’ command with the ‘DET’ mode option in FATCAT to view the quality of TORTOISE preprocessing and the gradient matching. For this, default tracking parameters were used (FA > 0.2, turning angle < 60°., and keeping tracts with length > 20 mm) to find tracts within the whole brain mask with OR logic. The resulting file was viewed in SUMA (Saad and Reynolds, 2012) along with a volume image to ensure that no obvious regions had been masked out and that major bundles appeared where expected. Uncertainty intervals (i.e., variances) of the FA and first eigenvector were estimated with ‘3dDWUncert’ using 500 jackknife-resampling iterations, for use in the probabilistic tractography.

Rs-fMRI data

The rs-fMRI data were analyzed using the Analysis of Functional Neuroimages (AFNI) software (Cox, 1996). The afni_proc.py script for preprocessing and @ANATICOR script for noise detection and removal from resting-state time series were used (Jo et al., 2013; Tinaz et al., 2016). This approach reduces the sensitivity of correlation coefficients to head motion and diminishes hardware-related noise more efficiently in resting-state data (Jo et al., 2010). The first two EPI volumes were removed to ensure that all remaining volumes were at magnetization steady-state. Spikes were identified in the time series. EPI volumes were slice-time corrected and motion parameters were estimated using rigid body transformations (three translations and three rotations, and their derivatives). The motion limit was 0.4 mm and outlier limit was 0.1. Subsequently, anatomical images were spatially normalized to a Montreal Neurological Institute (MNI) template (MNI_caez_N27). All transformations were applied at once to the EPI data to prevent multiple resampling steps. The EPI volumes were smoothed with a 6 mm full-width half-maximum Gaussian kernel. Nuisance variables (motion, spikes, local white matter) were regressed out. Global signal was not removed to prevent the introduction of spurious (anti)correlations (Saad et al., 2012). The time series were bandpass-filtered (0.01 < f < 0.1 Hz) to capture the resting-state fluctuations of the blood oxygenation level-dependent (BOLD) signal.

Data analyses

Network definition and metrics

We used the Craddock atlas (Craddock et al., 2012) with 200 cortical, subcortical, and cerebellar parcellations as our network template (Fig. 1). This atlas was created using rs-fMRI time series and yielded spatially coherent and temporally homogeneous clusters. We chose the option with 200 parcellations because this number was found to provide regional interpretability of the results without compromising cluster homogeneity or causing loss of substantial information. This template (from here on “Craddock200”) was applied to both DTI and rs-fMRI datasets, and every parcel in the template was considered a “node” in the network. The center-of-mass coordinates and labels of the nodes can be found in Supplementary material Table S1.

Fig. 1

Sample illustration of the Craddock parcels (n = 200) displayed on axial sections of the MNI_caez_27 template.

The characteristics of the graph metrics we used are defined as follows (Bullmore and Sporns, 2009; Rubinov and Sporns, 2010):

Node strength indicates how strongly one node is connected to the rest of the nodes in the network. We used the undirected, weighted sum of the connections.

Clustering coefficient measures the density of connections between neighboring nodes and is associated with the efficiency of local information transfer. Local efficiency is related to the clustering coefficient and reflects how relevant a node is for the communication among neighbors. Modularity is a measure of the community structure of a network. Communities are defined as groups of densely interconnected nodes that are sparsely connected with the rest of the network (Bertolero et al., 2015; Meunier et al., 2010; Mišić and Sporns, 2016; Newman, 2006). The modular organization of networks determines the network dynamics and information flow.

Tissue segmentation and ROI setup

The following steps were performed to bring the Craddock200-based parcellations in the MNI space into the subject-specific DTI space:

The T1-weighted volume was intensity-normalized using the 3dUnifize command in AFNI, and then skull-stripped in FreeSurfer (http://surfer.nmr.mgh.harvard.edu/). The skull-stripped T1-weighted volume was then co-registered to the MNI_caez_N27 template in AFNI using the auto_warp.py nonlinear registration function.

The skull-stripped T1-weighted volume was also aligned to the target T2-weighted volume using AFNI's align_epi_anat.py script with the local Pearson's coefficient (LPC) cost function (Saad et al., 2009).

This T1/T2 linear transformation was concatenated and applied to bring the Craddock200 template (originally in MNI space) to the DTI space using the 3dNwarpApply nonlinear transformation function.

In the end, an ROI volume with 200 parcellations was obtained for each subject which was subsequently used to perform ROI-based probabilistic tractography.

Probabilistic tractography

We used FATCAT for full probabilistic tractography. FATCAT uses repeated iterations of whole brain tracking to estimate the likelihood of structural white matter connections between all pairs of target ROIs within a network. Probabilistic tractography on the whole brain Craddock200 parcellation was performed using the command 3dTrackID with the PROB mode option. Tracking parameters used were: FA > 0.2, turning angle < 60°., keeping tracts with length > 20 mm, thresholding fraction > 0.021, five seeds per voxel and a total of 5000 Monte Carlo iterations. At each iteration, the voxel parameters were perturbed according to estimated parameter uncertainties. Tracks passing through individual ROIs and locations of tracts that intersected any pairs of ROIs were recorded. The probability statistics of the tractographic connections between all pairs of ROIs were automatically calculated. The automatically generated ‘*.grid’ file output of 3dTrackID contained matrices (here, 200 × 200) of output statistics for each subject describing the properties of the tractographic connections between all pairs of nodes, such as the number of tracts (NT) and the mean and standard deviation of FA. The probabilistic number of tracts between the nodes divided by the total number of probabilistic tracts in the whole brain was defined as the fractional NT (fNT) and used as the measure of anatomical pairwise connectivity between nodes. A 200 × 200 adjacency matrix for each subject was created based on the pairwise fNT values. As a confirmation, the white matter “degree” of each node was computed using the FA values and reported in Supplementary material.

White matter segmentation

We extracted the total white matter and white matter hypointensity volumes of each subject from the FreeSurfer segmentation files and normalized them to the total brain volume. The hypointensity volumes were used to identify the burden of nonspecific white matter disease.

The @ROI_Corr_Mat function in AFNI was employed to create the adjacency matrices for the subsequent network analyses. The Craddock200 template was used as the ROI volume from which the simple averages of the BOLD signal time courses were extracted for each node. A Pearson correlation coefficient was calculated between the average BOLD time course of each node and that of every other node creating a 200 × 200 adjacency matrix for each subject. The coefficients were then Fisher z-transformed.

Graph analyses

Graph metrics were calculated using custom MATLAB scripts containing functions from the Brain Connectivity Toolbox (BCT) (Rubinov and Sporns, 2010). The 200 × 200 adjacency matrices generated from DTI and rs-fMRI data were entered in graph analyses. For each node of the network, strength as a general measure of connectivity, and clustering coefficient and local efficiency as measures of local network properties were calculated. The DTI adjacency matrices were not thresholded because the fNT represents (the probability of) normalized white matter connectivity between nodes. For the rs-fMRI adjacency matrices, a range of thresholds was used (0.1 ≤ r ≤ 0.4 in 0.05 increments, total of seven).

Between-group differences in these graph metrics were assessed using permutation testing with 5000 iterations and a significance level set at p < 0.05, two-tailed (Mattest function in Bioinformatics toolbox of Matlab 2013a). Since threshold values are often arbitrarily determined (Rubinov and Sporns, 2010), for consistency, we chose stringent criteria for the rs-fMRI data and considered graph metrics that indicated a significant difference (p < 0.05) between the two groups at a minimum of five out of seven correlation thresholds.

Network modularity was calculated on group-averaged and weighted adjacency matrices using the multiscale community louvain algorithm in the BCT. This function provides optimal community structures by subdividing the networks into non-overlapping groups of nodes which maximizes the number of within-group edges, and minimizes the number of between-group edges. The resolution parameter gamma has a default value of 1, and values > 1 detect smaller modules. Here, a gamma range between 1 and 2 in 0.1 increments (500 iterations per increment) was used. The partition distance function was used to assess the similarity of modules across the gamma range. Higher similarity indicates stronger stability of the modules. In addition, for each gamma, groupings of at least three nodes were considered modules.

Statistics

Demographic data

Age, MMSE, and MoCA scores were compared between the groups using two-tailed, two-sample t-tests (p < 0.05). UPDRS part III (motor exam) “on” and “off” scores were compared using a two-tailed, paired-sample t-test (p < 0.05).

White matter segmentation

The normalized total white matter and white matter hypointensity volumes were compared between the two groups using two-tailed, two-sample t-tests (p < 0.05).

Post-hoc analyses

Structure — function overlap

We examined whether the functional reorganization in the PD connectome was supported by structural changes. As described in Section 2.5.2.5 in the functional graph analysis, first we identified the groups of nodes that showed significantly reduced values in each graph metric category (node strength, clustering coefficient, and local efficiency) in PDs compared to HVs. Then, the average nodal fNT values for each group separately was calculated per subject. Finally, we performed two-sample t-tests to compare these average nodal fNT values between the two groups (two-tailed, p < 0.05).

We also calculated the average fNT values of nodes belonging to each rs-fMRI module (at gamma = 1.3, 1.4, and 1.5. See in Results section) for each subject. Two-way repeated measures ANOVA test was performed using group and module type as the within-subject factors and average fNT value as the dependent variable.

Relationship between clinical and modularity data

To investigate the relationship between the clinical data and modular organization, we computed the average node fNT per DTI module and average unthresholded node strength per rs-fMRI module for each PD patient. These values corresponded to the average nodal white matter and functional connectivity strength of each module, respectively. We performed a separate multiple regression analysis for each of the DTI and rs-fMRI modularity datasets, in which the respective average fNT and unthresholded node strength values per module were entered as independent variables. We entered the UPDRS part III (motor exam) “on” scores as an objective measure of disease severity and the MoCA scores as a measure of cognitive function as the dependent variables in the analyses. We used the MoCA scores for two reasons: 1) most patients showed ceiling effect in the MMSE and 2) MoCA is more sensitive to executive dysfunction.

Results

Demographic and clinical data are summarized in Table 1. There were no significant differences in age (p = 0.8), MoCA (p = 0.5) and MMSE scores (p = 0.4), or gender and handedness between the PD and HV groups. There was 30% difference between the UPDRS part III (motor exam) “off” and “on” scores in the PD group which was also statistically significant (p < 0.0001). Since the “on” exam was performed immediately upon completion of the rs-fMRI scans, this difference further validates that the rs-fMRI scans were collected during the clinical “on” state.

Table 1

Demographic and clinical data.

	PD (n = 20)	HV (n = 20)
Age	62.5 ± 6.9	61.9 ± 6.6
Gender (M:F)	9:11	9:11
Handedness (L:R)	1:19	1:19
MoCA	27.6 ± 2.9	27.1 ± 2.1
MMSEa	29.4 ± 1.1	29.7 ± 0.6
Onset side (L:R)	9:11	–
Disease stage (H&Y)	2.0 ± 0.5	–
Disease duration	7.1 ± 3.3	–
UPDRS total
Off	45.4 ± 14.6	–
On	37.3 ± 12.6	–
UPDRS I	2.6 ± 1.7	–
UPDRS II	10.8 ± 4.6	–
UPDRS III		–
Off	29.5 ± 9.6	–
On	20.9 ± 7.9	–
UPDRS IV	3.0 ± 2.2	–

HV: healthy volunteer, H&Y: Hoehn & Yahr, MMSE: Mini Mental State Examination, MoCA: Montreal Cognitive Assessment test, PD: Parkinson's disease, UPDRS: Unified Parkinson's Disease Rating Scale.

The MMSE score of one HV subject is missing.

DTI data

Local network metrics

The DTI graph analysis revealed no significant difference in the average fNT values across the whole connectome between the two groups (p = 0.52). However, we observed reduced node strength, clustering coefficient, and local efficiency in a small group of nodes mostly in the frontoparietal regions and subcortically involving the caudate, hippocampus, and hypothalamus on the right in the PD group compared to HVs. Table 2 shows the MNI coordinates and labels of the nodes that showed significant differences in DTI graph metrics between the two groups.

Table 2

DTI graph metrics.

	x	y	z	NS	CC	LE
HV > PD
Parietal
R precuneus	1	− 55	19	–	✓	✓
L SMG	− 47	− 51	54	–	✓	✓
R post-CG	45	− 23	58	✓	–	–
	28	− 47	68	–	–	✓
Frontal
R SMA	5	− 1	72	–	–	✓
	1	21	63	✓	–	–
L MFG	− 40	24	45	–	✓	✓
	− 52	31	25	✓	–	–
L pre-CG	− 47	13	30	✓	–	–
Basal ganglia
R caudate	5	9	− 5	✓	–	–
Subcortical
R hypothalamus	1	− 1	− 6	✓	–	✓
R hippocampus	24	− 37	− 3	–	–	✓
PD > HV
None	–	–	–	–	–	–

MFG: middle frontal gyrus, SMA: supplementary motor area, SMG: supramarginal gyrus, pre/post-CG: pre/postcentral gyrus, L: left; R: right. NS: node strength; CC: clustering coefficient; LE: local efficiency.

Modularity

The overall degree of modularity (Q) was not significantly different between the groups across the whole gamma range for the DTI data (p = 1). Table S3a in Supplementary material shows the degree of similarity and numbers of the modules across all gamma thresholds for DTI data. Each cell in the matrix reflects the similarity value between modules. Both groups had very similar average similarity values across all thresholds. The highest average similarity value was 0.91 at gamma = 1.4, 1.5, and 1.6 for both groups.

The modular organization of the DTI connectome was virtually identical for both groups across all gamma thresholds. At gamma = 1.5, it consisted of eight modules for both groups: module 1: left hemisphere occipito-temporal and basal ganglia, module 2: medial frontoparietal, module 3: left hemisphere orbito-temporal, module 4: cerebellar, module 5: right hemisphere occipito-temporal, module 6: right hemisphere basal ganglia, module 7: left hemipsphere temporo-parieto-frontal, module 8: right hemisphere frontoparietal (Fig. 2). Tables S4a and b in Supplementary material lists the node compositions of the DTI modules.

Fig. 2

The DTI modules for both groups at gamma = 1.5 are shown as circular graphs using Gephi (Jacomy et al., 2014). In both groups, the same arbitrary threshold for white matter connection strengths was used for display purposes. BG: basal ganglia, Cb: cerebellum, Fron: frontal, Occ: occipital, Orb: orbital, Par: parietal, Temp: temporal. L: left, R: right, m: medial. The node contents of modules are listed in Table S4b in Supplementary material.

White matter segmentation

There were no significant differences between the two groups in normalized total white matter and white matter hypointensity volumes. The average total white matter volume was 516,250.9 ± 75,693.3 mm3 in the PD and 512,822.1 ± 70,813.4 mm3 in the HV group (p = 0.67). The average white matter hypointensity volume was 3301.7 ± 2651.6 mm3 in the PD and 3385.9 ± 4972.5 mm3 in the HV group (p = 0.90) (also see Supplementary material Table S2).

Rs-fMRI data

The rs-fMRI graph analysis also revealed no significant difference in the average unthresholded node strength of the whole connectome between the two groups (p = 0.15). However, we found significantly reduced node strength in 31 nodes, clustering coefficient in 55 nodes, and local efficiency in 52 nodes diffusely distributed across the whole connectome including the striatal nodes bilaterally in the PD group compared to HVs. Table 3 shows the MNI coordinates and labels of the nodes that showed significant differences in rs-fMRI graph metrics between the two groups.

Table 3

Rs-fMRI graph metrics.

	x	y	z	NS	CC	LE
HV > PD
Occipital
L MOccG	− 50	− 84	− 1	✓	✓	✓
Temporal
L ITG	− 45	− 65	− 8	✓	✓	✓
R MTG	63	− 51	16	✓	–	✓
	59	− 27	− 8	✓	✓	✓
	63	− 49	− 1	–	✓	✓
	63	− 20	− 16	–	✓	✓
L MTG	− 60	− 48	5	✓	✓	✓
	− 60	− 31	− 1	–	✓	–
	− 60	− 12	− 23	–	✓	✓
	− 60	− 48	− 6	–	✓	✓
	− 60	4	− 32	–	✓	–
	− 60	− 68	16	–	✓	✓
R STG	64	− 12	5	–	✓	✓
L STG	− 60	− 26	8	–	–	✓
R TP	60	8	− 8	✓	✓	✓
L TP	− 49	9	− 8	✓	✓	✓
	− 27	10	− 23	–	✓	✓
R PHG/amy	24	2	− 25	–	✓	✓
R fusiform	43	− 13	− 32	–	✓	–
Parietal
R AG	32	− 66	49	✓	✓	✓
L AG	− 48	− 65	49	✓	✓	✓
R SPL	24	− 63	73	–	✓	✓
	47	− 45	60	–	✓	✓
L SPL	− 28	− 74	55	✓	✓	✓
R SMG	63	− 33	27	✓	✓	✓
	59	− 33	51	–	✓	–
	49	− 30	27	✓	✓	✓
	49	− 69	49	✓	✓	✓
	60	− 54	41	–	✓	–
L SMG	− 60	− 48	31	–	✓	✓
	− 47	− 51	54	–	✓	–
	− 60	− 25	39	–	✓	–
L post-CG	− 62	− 9	25	–	✓	✓
Frontal
R MOrbG	31	57	− 17	✓	✓	✓
	31	47	− 22	✓	✓	✓
L MOrbG	− 50	54	− 5	✓	✓	✓
	− 27	57	− 17	–	✓	✓
R SMA	6	0	57	–	✓	✓
	1	20	55	✓	✓	✓
R SFG	47	27	25	✓	✓	✓
R MFG	44	50	1	✓	✓	✓
	47	26	44	–	✓	–
L MFG	− 28	18	59	✓	✓	✓
	− 50	42	9	✓	✓	✓
R IFG (p. Op)	53	13	27	✓	✓	✓
L IFG (p. Op)	− 50	10	19	✓	✓	✓
R IFG (p. Orb)	37	27	− 8	✓	✓	✓
	28	28	− 22	✓	✓	✓
	58	24	− 8	✓	✓	✓
	47	29	− 14	–	✓	✓
L IFG (p. Orb)	− 49	34	− 8	✓	✓	✓
R ACC	4	14	25	–	✓	–
L ACC	− 3	37	25	–	✓	✓
L OFI	− 33	24	− 8	✓	✓	✓
L post cing	− 0	− 19	36	✓	✓	✓
L pre-CG	− 60	12	34	–	✓	✓
Basal ganglia
L caudate	− 17	11	15	✓	–	✓
R caudate	13	10	6	✓	✓	✓
	5	9	− 5	–	–	✓
R putamen	26	7	− 2	–	–	✓
PD > HV
Cerebellum
L Cb (culmen)	− 28	− 42	− 31	✓	–	–
L Cb (vermis)	− 0	− 35	− 24	✓	–	–

ACC: anterior cingulate cortex; AG: angular gyrus; Amy: amygdala; Cb: cerebellum; Cing: cingulate; IFG, MFG, SFG: inferior, middle, superior frontal gyrus (p: pars, Op: opercularis; Orb: orbitalis); OFI: orbitofrontal insula, SPL: superior parietal lobule; ITG, MTG, STG: inferior, middle, superior temporal gyrus; MOccG: middle occipital gyrus; MOrbG: middle orbital gyrus; PHG: parahippocampal gyrus; pre/post-CG: pre/postcentral gyrus; SMA: supplementary motor area; SMG: supramarginal gyrus; TP: temporal pole. L: left; R: right; Post: posterior. NS: node strength; CC: clustering coefficient; LE: local efficiency.

The overall degree of modularity (Q) was also not significantly different between the groups across the whole gamma range for the rs-fMRI data (p = 0.46). Table S3b in Supplementary material shows the degree of similarity and numbers of the modules across all gamma thresholds for rs-fMRI data. Each cell in the matrix reflects the similarity value between modules. Both groups had very similar average similarity values across all thresholds. The highest average similarity value was 0.91 at gamma = 1.4, 1.5, and 1.6 for both groups. The HV group had the highest average similarity value (0.76) for modules at gamma = 1.5, whereas the PD group the highest average similarity value (0.82) for modules at gamma = 1.4.

Both groups showed a similar modular organization at gamma = 1.3, and the following modules were identified: module 1: parietal and premotor nodes (sensorimotor network); module 2: temporal and ventral prefrontal nodes (task-set maintenance network); module 3: medial and lateral temporal and basal ganglia nodes (temporal network); module 4: fronto-parietal nodes (dorsal attention network); module 5: occipital nodes (visual network); module 6: cerebellar and basal ganglia nodes; module 7: medial fronto-parietal nodes (default mode network). The HV group had an additional module of orbitofrontal and temporal nodes at this threshold.

At higher thresholds, the HV group demonstrated a more fine-grained modular organization compared to the PD group. For instance, the HV group continued forming more submodules at gamma = 1.5: The sensorimotor network was further divided into medial and lateral parietal/premotor modules. The dorsal attention network was divided into premotor/prefrontal and lateral frontoparietal modules. The default mode network was divided into medial frontoparietal and dorsomedial prefrontal/temporal modules. There were also additional orbital and basal ganglia modules. On the other hand, PD group showed two additional small modules at this threshold: lateral frontal and medial parietal (Fig. 3).

Fig. 3

The rs-fMRI modules for both groups at gamma = 1.5 are shown as circular graphs using Gephi (Jacomy et al., 2014). In both groups, the same arbitrary threshold for functional connection strengths was used for display purposes. BG: basal ganglia, Cb: cerebellum, Fron: frontal, Occip: occipital, Orb: orbital, Par: parietal, PF: prefrontal, PM: premotor, Temp: temporal, d: dorsal, l: lateral, m: medial, v: ventral. The node contents of modules are listed in Table S5c in Supplementary material.

Tables S5a and b in Supplementary material list the node compositions of the rs-fMRI modules for both groups at gamma thresholds of 1.3 and 1.5, respectively. The nodes that were shared by both groups in each module and the different ones are listed separately.

Post-hoc results

Structure — function overlap

There were no differences in the average fNT values between the two groups in the nodes that showed significantly reduced functional node strength (p = 0.92), clustering coefficient (p = 0.74), and local efficiency (p = 0.66) in the PD group compared to HVs. In other words, the nodal functional connectivity differences between the groups were not supported by structural connectivity differences.

The two-way repeated measures ANOVA comparing the average fNT values of rs-fMRI modules between the groups did not reveal any main effect of or interaction between group and module type at any of the gamma thresholds 1.3, 1.4, or 1.5.

Relationship between clinical and modularity data

We did not observe any correlations between the clinical data and modularity findings. Motor severity and cognitive scores (UPDRS part III “on” and MoCA) did not correlate with the modular node strength in the rs-fMRI data (tested at gamma thresholds 1.3 and 1.5) or in the DTI data (tested at gamma = 1.5).

Discussion

The connectome, both in its structure and function, is altered in non-demented PD patients on medication. Here, we summarize and discuss the results.

White matter connectivity

The PD group showed significantly decreased white matter connectivity in the right caudate and hypothalamus, and in a number of frontoparietal nodes, compared to HVs. Most of these frontoparietal nodes and the right hypothalamus and hippocampus also showed reduced local connectivity.

White matter lesions of ischemic etiology are commonly observed in clinical scans of older individuals and are not routinely reported or accounted for in DTI analyses. Our cohorts were well-matched in global white matter and white matter lesion volumes. Therefore, we think that the differences in white matter connectivity between the groups cannot be explained based on nonspecific white matter disease. Moreover, the anatomical pattern of reduced white matter connectivity is also consistent with disease-specific pathology in PD. Our findings are in line with reports of impaired white matter integrity in animal models and in vivo imaging studies in PD. The DTI studies to date demonstrated significantly reduced FA values in major tracts (e.g., parietal and frontal white matter) in PD patients with minimal cognitive impairment or dementia, but not in those with intact cognition (Hattori et al., 2012; Melzer et al., 2013). Of note, our PD cohort was not demented and was well-matched to the HV group in cognitive performance. The DTI and graph theory-based methods used here enabled us to detect local changes in white matter integrity even at relatively early stages of the disease process in PD.

Contrary to our hypothesis, the PD group showed virtually the same modular organization of white matter connectivity as the HV group which remained stable across multiple resolution thresholds. This indicates that despite subtle connectivity changes, the overall structural organization of the PD connectome remains robust. However, in later stages of the disease process more pronounced loss of white matter connections might occur and lead to changes in the modular architecture. This possibility needs to be investigated in longitudinal studies.

Functional connectivity and the probable role of dopaminergic treatment

The PD group showed significantly decreased node strength compared to HVs in 31 nodes distributed across the cortex. About twice as many nodes also showed significantly reduced local connectivity. Critically, striatal nodes also exhibited these connectivity changes in the PD group. The PD group also showed reduced functional modularity compared to HVs across higher resolution thresholds. More specifically, the modular delineation was reduced in the PD group compared to HVs in the sensorimotor and two core cognitive networks, namely the dorsal attention and default mode networks, as well as in the orbitofrontal cortex and basal ganglia.

Several rs-fMRI studies in PD reported reduced functional connectivity in various networks. Some studies also demonstrated correlations between the connectivity findings and disease severity and cognitive status. There is considerable variability in the results of these studies that seems to stem from the differences in analysis methods, cohort characteristics, disease severity, cognitive status (Baggio et al., 2014, Baggio et al., 2015, Lebedev et al., 2014), and medication state. That said, there is convergence between our findings and those of other studies that used similar cohorts. For instance, in a study with non-demented PD patients on medication compared to matched controls, an aberrant positive resting-state correlation between the right central executive and default mode networks was found (Putcha et al., 2015). We think that the relatively reduced modularity in the dorsal attention and default mode networks in our PD group is in line with this finding. By definition, reduced modularity of a network implicates aberrant connections with other networks.

We think that medication status during scanning also plays a crucial role. In fact, a recent study investigated the functional reorganization of ten brain networks composed of 226 nodes in 19 cognitively intact PD patients on and off medication (Berman et al., 2016). PD patients off medication did not show a significant change in global efficiency or overall local efficiency compared to controls. On the other hand, l-dopa intake led to overall reduced local efficiency, most significantly in the dorsal attention network in PD patients (PD-on compared to PD-off). Compared to controls, PD patients on medication showed a significant decrease in local efficiency in the sensorimotor and a significant increase in local efficiency in the subcortical network (Berman et al., 2016). The PD cohort in the study by Berman et al. was similar to ours in size and disease characteristics, but the graph analysis method was different (unweighted binarized graphs in Berman et al. as opposed to weighted graphs here) which could explain the discrepancies in our findings.

Nevertheless, it is important to note that despite methodological differences between studies in PD, the results indicate that dopaminergic treatment does not improve or normalize the functional network properties unanimously to the degree of those observed in matched controls (Berman et al., 2016, Göttlich et al., 2013, Putcha et al., 2015, Tessitore et al., 2012). We think that our rs-fMRI nodal and modular findings are also related to these differential effects of dopamine.

It is known that dopamine exerts a dose-dependent influence on performance in different tasks in PD. For example, performance in stimulus-reward association learning tasks that recruit the ventral striatum and orbitofrontal cortex worsens with dopaminergic treatment in PD (Cools et al., 2002, Cools et al., 2006, Peterson et al., 2009). On the other hand, performance in cognitive tasks that recruit the dorsal striatum and dorsal frontoparietal cortex (e.g., working memory, planning, task-switching), improve with dopaminergic treatment in PD (Cools et al., 2001, Macdonald and Monchi, 2011). According to the “dopamine overdose” hypothesis, this discrepancy can be explained by the different degrees of pathology affecting the dorsal and ventral parts of the striatum and the associated cortical circuits (Kish et al., 1988). The more affected dorsal part is more severely deprived of dopamine and functions more optimally with dopamine replacement, whereas the less severely affected ventral part is overwhelmed by dopamine replacement and functions suboptimally (Cools et al., 2001, Cools et al., 2003, Gotham et al., 1988). However, the dose-dependent influence of dopamine is not necessarily linear and dopamine levels can also exhibit an inverted U-shaped relationship between behavioral performance and BOLD response (Cools, 2006, Rowe et al., 2008).

Neuroimaging studies have also provided evidence for differential effects of dopaminergic states on striatal functional connectivity patterns. In a study with drug naive de novo PD patients, the striatal dopamine levels correlated differentially with the resting-state whole-brain functional connectivity patterns of the caudate and putamen seeds. For instance, functional connectivity between the anterior putamen and mesial frontal areas correlated negatively with striatal dopamine levels, whereas functional connectivity between the dorsolateral prefrontal areas and the anterior and posterior putamen correlated positively with striatal dopamine levels (Baik et al., 2014).

Taken together, dopaminergic treatment might differentially influence various frontoparietal-striatal circuits. It is conceivable that this differential influence results in less correlated BOLD signal between these regions and the rest of the network as our reduced nodal functional connectivity findings indicate. Moreover, reduced modularity in the dorsal attention, sensorimotor, and default mode networks as well as lack of modular distinction in the orbitofrontal and basal ganglia nodes in the PD group also implicate loss of local specificity in the community organization of these networks probably due to less correlated BOLD signal among neighboring nodes. Nevertheless, the differential influence of dopamine is likely determined by additional factors. Inter-individual differences and disease severity would also be expected to effect the neurovascular response to dopamine. By itself, this study cannot determine whether these connectivity differences are due to the disease state itself or the interaction between medication and the disease state. Therefore, the specific role of dopaminergic treatment in functional modularity should be tested directly by assessing the same PD cohort on- and off-medication. Moreover, in addition to the neuronal effects, the direct vascular effects of dopamine in neurovascular coupling might also be an important factor (Choi et al., 2006; Zaldivar et al., 2014) affecting the functional connectivity patterns.

Structure-function relationship

There is evidence from experimental and network modeling studies suggesting that structural connectivity shapes and constrains functional connectivity (Hagmann et al., 2010, Hermundstad et al., 2013, Honey et al., 2007, Honey et al., 2009). The results of our study and others add another layer of complexity to this organizational principle. We found many more nodes that showed significant functional as opposed to structural connectivity differences in the PD group in the “on” state, and they were not precisely co-localized. In addition, at the large-scale network level, the functional modularity changes in the PD group were not supported by structural modularity. Another study in patients with schizophrenia demonstrated altered functional connectivity in three spatially distinct large-scale networks compared to controls (Cocchi et al., 2014). Co-localized alteration of functional and structural connectivity was found in two of the three networks, yet the relationship was in opposite directions, namely, reduced structural and functional connectivity in one network, and reduced structural but enhanced functional connectivity in the other (Cocchi et al., 2014).

These observations indicate that the structure-function relationship in the connectome of patient populations may be nonlinear. A possible explanation for this nonlinear relationship specifically in PD could be synaptic dysfunction and changes in synaptic plasticity. Growing evidence suggests that neurodegeneration in PD is more likely to start at the synapse (Bellucci et al., 2016, Picconi et al., 2012; Schulz-Schaeffer, 2010). The highest percentage of alpha-synuclein aggregates in PD are localized at the presynapses. Several other PD-associated proteins (e.g., LRRK2, parkin, DJ-1, PINK1) also alter the dopaminergic presynaptic site (Beccano-Kelly et al., 2015; Kitada et al., 2009; Madeo et al., 2014; Scherfler et al., 2004). The presynaptic changes and depletion of dopamine in turn alter the postsynaptic plasticity (Day et al., 2006). These changes lead to severe synaptic dysfunction that may result in retrograde axonal damage and ultimately neuronal cell body death (Anichtchik et al., 2013, Zaltieri et al., 2015). It is plausible that in the relatively early disease stages, synaptic dysfunction results in relatively subtle axonal loss without a significant change in structural connectivity pattern, but leads to more robust functional connectivity changes by altering the BOLD signal as we observed in this study. This hypothesis needs to be tested further in longitudinal studies.

Relationship between modularity and behavior

Finally, in our post-hoc analyses, we did not observe any correlation between the MoCA and UPDRS part III (motor exam) scores and rs-fMRI and DTI modularity findings in PD. Since the DTI modularity was unaffected in PD, the lack of a structure-behavior correlation is probably not surprising. As for the lack of a function-behavior relationship; the functional modularity probably reflects a complex combination of disease-related long-term and medication-related transient (re)organization, and may thus not correlate with rather temporary measures of cognitive and motor function. We think that the cognitive and motor consequences of the functional modularity changes warrant further investigation in larger patient cohorts on- and off-medication, stratified based on the level of motor dysfunction and cognitive deficits.

In conclusion, our results provide additional evidence regarding changes in white matter and functional connectivity in non-demented PD patients compared to matched healthy controls. The regional reduction in white matter connectivity, along with the reduction in regional functional connectivity and functional modularity in cognitive and motor networks, provide important information regarding the pathophysiology of PD. Our results can also be extended to monitoring the effects of therapeutic interventions on the whole connectome in PD.

Funding

This work was supported by the National Institute of Neurological Disorders and Stroke Intramural Research Program.