Network structure of brain atrophy in de novo Parkinson's disease

1Parkinson's disease (PD) is the second most common neurodegenerative illness. Its clinical hallmarks are due to loss of dopamine neurons in the substantia nigra (SN); however, post-mortem studies have shown that PD pathology affects several other brain areas. Surprisingly, however, magnetic resonance imaging (MRI) studies have failed to consistently demonstrate regional brain atrophy at least in the earlier stages of the disease. This is possibly due to studies consisting mostly of small numbers of subjects from single centers and the use of methodological tools that are relatively insensitive to subcortical atrophy, coupled with univariate methods that suffer from lack of statistical power.

2The distribution of Lewy body pathology in post-mortem samples led Braak et al. to propose that a disease process spreads from the brainstem to subcortical areas and then to the cerebral hemispheres along neuronal pathways. More generally, the network-spread hypothesis suggests that all neurodegenerative diseases target spatially distributed intrinsic brain networks. At a molecular level, neurodegenerative diseases are now thought to involve prion-like spread of toxic misfolded protein aggregates. Alpha-synuclein fibrils, implicated in PD pathogenesis, have recently been shown to spread from cell to cell via neuronal pathways after inoculation in mouse brain.

3Neuroimaging studies have revealed that the spatiotemporal organization of the brain takes the form of intrinsic connectivity networks (ICNs), which are consistent across studies, in test–retest analysis, and during both rest and task states. Intrinsic networks have also been thought to be potential targets for AD since the earliest reports of default mode network dysfunction in the disease, and analysis of large MRI data sets in AD and other dementias have supported the network-vulnerability hypothesis; however, it has yet to be tested in PD. Here, we analyzed data from newly diagnosed PD patients (n = 232) and an age-matched control group (n = 117) obtained from the Parkinson's Progression Markers Initiative (PPMI) database (www.ppmi-info.org/data), an observational, multicenter longitudinal study designed to identify PD progression biomarkers. We used deformation-based morphometry (DBM) and tensor probabilistic independent component analysis (ICA) to identify brain regions demonstrating atrophy in early PD.

4We also sought to provide support for the network-spread hypothesis in PD by showing that the disease, in humans, targets intrinsic brain networks. We compared the set of atrophic regions in PD patients to ICNs from young healthy subjects and tested them for spatial overlap. To further investigate the model of disease propagation through brain networks, we used resting state functional MRI (fMRI) and diffusion-weighted MRI (DW-MRI) of healthy subjects to define the normal brain connectome and determined whether the pattern of atrophy in PD was compatible with spread via this network from a presumed disease epicenter in the SN.

5DBM was used as the measure of local brain atrophy. It is a measure of the change in the shape of each voxel that results from applying non-linear spatial normalization of the entire brain to a standard template. For each subject, we obtained one parametric image of MRI-determined voxel-wise difference in volume, compared to the template brain. ICA was conducted on these DBM images using FSL MELODIC software. No constraint was imposed on the number of components, and probabilistic ICA estimated 30 independent components of deformation in the PPMI data set (PD patients and controls combined). Each ICA component consists of a spatial map and the average deformation value in that map for each subject.

6In each of these 30 independent components, the average deformation between PD and control groups was compared using an unpaired t-test. PD patients had significantly lower DBM values in one and only one of the 30 deformation components (p = 0.0001 uncorrected, p = 0.003 with Bonferroni correction). The next spatial component in terms of statistical significance (p = 0.06, uncorrected) consisted of cerebellar areas VIIIA, VIIB, CrusII, and IX known for their involvement in motor and executive function. None of the other components demonstrated a difference between PD and controls (p > 0.05 uncorrected for all remaining components). We will call the deformation pattern showing a group difference the PD-ICA network from this point on. Ten of the other ICA components demonstrated an effect of age in the entire group. In three of these, there was a positive correlation between the component expression and age (meaning increased volume with age) and all three consisted of ventricle or cerebrospinal fluid space enlargement (Figure 1—figure supplement 1). The other seven age-related components demonstrated a negative correlation (volume loss with age) and consisted of areas of gray or white matter (Figure 1—figure supplement 2). The PD-ICA network also demonstrated greater atrophy with increasing age in both PD patients (r = −0.38, p < 10−9) and controls (r = −0.53, p < 10−9).

7Each spatial map was converted to a z-statistic image via a normalized mixture model and then thresholded at z ≥ 3. Regions were identified using the atlases of and. The PD-ICA network includes all components of the basal ganglia (substantia nigra, subthalamic nucleus, nucleus accumbens, putamen, caudate nucleus, and internal and external globus pallidus), the pedunculopontine nucleus, basal forebrain, including bed nucleus of the stria terminalis and an area containing the nucleus basalis of Meynert, the hypothalamus, amygdala, hippocampus, parahippocampal gyrus, and two thalamic regions, the ventrolateral nucleus and pulvinar. Cortical regions in this network are the insula, occipital cortex Brodmann area 19, superior temporal gyrus, rostral anterior cingulate cortex, premotor and supplementary areas, and parts of lateral prefrontal cortex (Figure 1, Tables 1, 2; see also Figure 1—figure supplement 4).

8To confirm that the PD-ICA network identified above was disease-related, we compared individual deformation values in the network to measures of disease severity. The two clinical measures used were the striatum binding ratio (SBR) measured with single photon emission computed tomography (SPECT) using the tracer [123I]FP-CIT to measure dopamine nerve terminal density, and the score on the Movement Disorder Society revised Unified Parkinson's Disease Rating Scale (UPDRS) part III, an objective measure of motor disability. For SBR, we used the average value of left and right putamen. There was a significant correlation between individual SBR and DBM values in the PD-ICA network in the PD group (r = 0.23, p < 0.0005, Figure 2). This shows that the greater the loss of dopamine nerve terminals, the greater the volume loss in the PD-ICA network. There was also a significant correlation between these two measures in the control group (r = 0.33, p < 0.0005).

9There was a significant correlation between DBM values within the PD-ICA network and UPDRS III in the PD patients (r = −0.22, p < 0.001; Figure 2). SBR was not significantly correlated with age in the PD subjects (r = −0.10, p = 0.12) but it was in controls (r = −0.35, p < 0.0001). Also, in PD subjects, SBR was significantly correlated with UPDRS III (r = 0.20, p = 0.002). We also tested whether UPDRS III was correlated with DBM values obtained from any one of the other 29 ICA components. There was only one other component marginally correlated with disease severity (r = −0.2, p = 0.048, Bonferroni corrected) consisting of the previously mentioned cerebellar network, areas VIIIA, VIIIB, CrusII, and IX. Because age, UPDRS III, and SBR all correlated with PD-ICA DBM, we performed multiple linear regressions. In the PD group, the model (PD-ICA ~ 1 + Age + UPDRS III + SBR) showed an effect of age (p = 2.4 e−08), UPDRS III (p = 0.06) and SBR (p = 0.01). In the controls, the model (PD-ICA ~ 1 + Age + SBR) demonstrated an effect of age (p = 4.7 e−08) but not SBR (p = 0.17).

10Finally, multivariate analysis was used to look for an effect of scanning site. For each obtained DBM-network, we applied the model: DBM ~ Group (PD/Control) + Age + Gender + Site. There was no significant effect of site after correcting for multiple comparisons (p > 0.1).

11We next tested the hypothesis that the PD-ICA deformation network represents an intrinsic functional network. We compared the PD-ICA network as well as the other 29 ICA maps obtained from the DBM ICA analysis to intrinsic brain networks in healthy brain. In order to increase confidence in the results, we determined normal intrinsic brain networks in several different ways:A seed-based functional connectivity map obtained from resting-state fMRI (rsfMRI) in 51 healthy volunteers with an a priori region of interest in the SN.Two sets of 70 and 100 typical functional resting-state networks from healthy volunteers previously identified with MELODIC ICA by and.

12A seed-based functional connectivity map obtained from resting-state fMRI (rsfMRI) in 51 healthy volunteers with an a priori region of interest in the SN.

13Two sets of 70 and 100 typical functional resting-state networks from healthy volunteers previously identified with MELODIC ICA by and.

14First, we compared the resting-state SN seed-based map (Figure 3B) with all 30 structural maps obtained from the PPMI ICA analysis. Only two ICA networks passed the threshold of |r| > 0.25. These two networks correspond to very similar subcortical networks in basal ganglia and brainstem, one of which is mostly white matter areas (r = 0.36), while the other is mostly gray matter areas (r = 0.30). The latter is the aforementioned PD-ICA network in which atrophy correlated with disease severity (Figure 3—figure supplement 1).

15Then, we compared the PD-ICA network to the 70 rsfMRI ICNs in normal brain provided by. One of the 70 networks passed the threshold (with r = 0.32). This network has been related to reward tasks, interoceptive functions, and motor/sensory processing (Figure 3C). We assessed statistical significance by generating 1000 permutations of each of the 70 ICNs by reassigning the coordinates of each voxel randomly (Figure 3—figure supplement 2). We then repeated this comparison using a finer decomposition of 100 resting-state ICNs from the Human Connectome Project (HCP) using MELODIC. Four components showed spatial overlap to the PD-ICA network using spatial correlation (Figure 3—figure supplement 3). The mean fMRI time series from these components were then used to determine whether they themselves belonged to one larger ICN. These time series demonstrated significant inter-correlation (p < 0.0016 by permutation testing). Finally, hierarchical clustering confirmed that all four components clustered together.

16We also compared the PD-ICA network to a map of regions responding to stimulus value during fMRI experiments as identified by meta-analysis. In these experiments, subjects typically evaluated an offered item (say a food) and experimenters identified brain regions where the fMRI signal tracked subjective value (e.g., willingness to pay for the item). The premise is that PD may affect dopamine projection sites that encode aspects of stimulus value. Only one of the 30 networks passed the threshold (r = 0.28), namely the aforementioned PD-ICA network (Figure 3D).

17In summary, the PD-ICA network exhibited significant spatial overlap with presumed intrinsic brain networks determined by three different methods.

18The sequential propagation model predicts that the spatial progression of the disease process will be determined by brain network topology. Connectivity between any region and the presumed disease epicenter will determine how severely it is affected. Here, we evaluated this assumption by exploring whether the gray matter atrophy patterns observed in PD patients could be explained by functional and geodesic distance (i.e., the number of edges separating two nodes in a graph) to the hypothetical pathogenic epicenter. We chose to use the SN as the epicenter based on known PD pathology. Note that the SN is unlikely to be the first affected site in the central nervous system (CNS); however, we postulate that it is likely to function as a source for propagation to the supratentorial brain. Network connectivity in health can be defined functionally, using rsfMRI, or structurally, using DW-MRI. The influence of disease on each node can be estimated from the statistical difference in deformation between PD and control groups. We parcellated the brain into 112 regions of interest (ROIs, Figure 4—figure supplement 1) and computed the degree of deformation in each region. We also generated two connectomes from these ROI using rsfMRI and DW-MRI data from two different pools of healthy subjects. The connection strength between each pair of regions was computed as described in the ‘Materials and methods’.

19There was a significant correlation between resting-state functional connectivity of each node with SN (in healthy brain) and the PD-related deformation (PD minus control t-score), (r = 0.40, p < 0.0001). We repeated the same analysis controlling for spatial proximity between each region and SN by entering Euclidean distance as a covariate. The correlation was unchanged (r = 0.38, p < 0.0001), suggesting that the relationship cannot be explained by spatial proximity. The correlation implies that higher functional connection between a given region and SN is related to higher PD-related atrophy in that region (Figure 4). When using an anatomical measure of connectivity (DW-MRI) (Figure 4), we also observed a significant relationship between the level of atrophy of each region and its geodesic distance to the SN (r = −0.28, p < 0.003). These results were not different after controlling for Euclidean distance to SN (r = −0.25, p < 0.005). Finally, we repeated this analysis using every ROI as a potential disease propagator and found that SN was the likeliest disease propagator when using the rsfMRI connectome (Table 3). However, the red nucleus (r = 0.28) and subthalamic nucleus (r = 0.28) were also identified as potential propagators. Repeating this analysis using a tractography-derived connectome also revealed that the SN was one of the likeliest propagators, but numerous cerebellar regions also emerged as potential propagators (Figure 4—source data 2). This may be due to difficulty in accurate identification of the targets of brainstem white matter tracts using DW-MRI.

20The combination of DBM and ICA of an unprecedentedly large set of MRI data at a magnetic field strength of 3T allowed us to map out the brain regions affected in de novo PD. Patients included in this study were diagnosed on average 7 months prior to the investigation and had no evidence of dementia (Table 4). Most MRI studies to date using T1-weighted images in PD have reported normal volumes of brainstem, basal ganglia, and cerebral cortex. In the SN, measurements of volume and shape have been inconclusive, reporting either no change, a decrease, or increase in volume depending on the method used. However, most studies using parametric mapping have failed to report a difference in SN. Reduced putamen or caudate volume has been reported but only in advanced cases with mild cognitive impairment or dementia. Similarly, hippocampus and amygdala atrophy are typically linked to cognitive impairment (for review see; ). With regard to cortical areas, studies using voxel-based morphometry (VBM), DBM, or cortical thickness measurements have reported significant differences between PD and controls in parahippocampal gyrus, inferior, middle and frontal temporal gyrus, parietal lobe and occipital lobe, but, once again, typically in advanced patients with cognitive impairment. Our finding of an atrophy network including brainstem and subcortical regions in these de novo patients is likely due to the larger number of participants available through the PPMI database compared to previous relatively underpowered investigations.

21Eidelberg and colleagues were the first to use data reduction techniques (principal component analysis) to map PD pathology using [18F]fluorodeoxyglucose (FDG) positron emission tomography. They identified a PD-related pattern (PDRP) implicating several regions identified in our analysis, namely globus pallidus, putamen, thalamus, premotor and supplementary motor areas. The expression of this PDRP differentiates patients from healthy controls and correlates with measures of disease severity. However, it probably reflects functional effects of the disease more than neuronal atrophy, as evidenced by the fact that it is normalized by levodopa or deep brain stimulation.

22We found that the set of regions demonstrating atrophy in the PD group corresponded spatially to one ICN in normal brain. rsfMRI and anatomical imaging have identified consistent sets of regions that act as functional networks, by virtue of anatomical connectivity and temporal covariance of neuronal activity. We showed that our PD-ICA set of regions corresponded to a normal ICN identified by seed-based resting-state functional connectivity and ICA. The finding that PD targets a set of connected brain regions supports the network-spread hypothesis, although an alternate explanation is that neurons in one network could share a common vulnerability. The brain regions that make up the PD-ICA network are involved in reward and motivation, as demonstrated by the fact that they respond to the subjective value of a perceived stimulus during fMRI. All the regions of this value network receive projections from midbrain dopamine neurons, in keeping with this neurotransmitter's role in signaling incentive value.

23Recent evidence supports a hypothesis originally put forward by Braak, according to which pathogenic forms of the protein alpha-synuclein spread throughout the nervous system leading to a stereotypical step-wise pattern of neurological impairment in PD. Misfolded synuclein fibrils injected focally into rodent brains spread to neuronally connected but not adjacent areas (see also Figure 1—figure supplement 3).

24Here, we tried to identify a putative best propagator or epicenter of the disease by creating functional and anatomical connectomes from an anatomical brain atlas. We showed that the SN could act as a propagator, since the rate of atrophy in any brain region was proportional to its effective topological distance from the SN, as determined by either functional (rsfMRI) or anatomical (DW-MRI) connectivity. We note that, while the SN contains dopamine neurons whose degeneration accounts for all the key clinical features of early PD, the Braak scheme identifies the dorsal motor nucleus (DMN) of the vagus in the medulla as the first affected CNS structure. Indeed, recent evidence supports the potential transfer of pathogenic alpha-synuclein from the intestine to the DMN via the vagus nerve. In the current study, the DMN did not demonstrate atrophy in PD compared to control subjects. Perhaps this is due to insensitivity of T1-weighted MRI-based anatomical methods in identifying neuronal loss in lower brainstem structures. The DMN is much smaller that the SN (approximately 1/10–1/20 in volume) and has limited contrast in T1 MRI. Our results, however, suggest that the SN could act as the propagator of the disease to the supratentorial CNS. Note, however, that spatial resolution limitations for all of the imaging modalities used here, T1-weighted MRI, fMRI, and DW-MRI, make it difficult to localize atrophy to small nuclei with complete certainty. Spatial inaccuracy in atlas generation and normalization of subject images to Montreal Neurological Institute (MNI) space may compound this problem. This is especially relevant to small structures of the basal forebrain and midbrain. For example, the propagator analysis identified SN, subthalamic nucleus, and red nucleus as potential propagators, but it is clear that the imaging techniques used here do not allow us to fully resolve these structures, either anatomically or functionally. Furthermore, because T1-weighted MRI is sensitive to iron content, changes in iron accumulation, for example, in SN, globus pallidus, or red nucleus, may be interpreted as volume changes by the DBM methodology.

25Nonetheless, comparing the atrophy pattern identified here to the stages described by,, we note that all the areas identified in stage 3 (pedunculopontine nucleus, amygdala, basal forebrain, and substantia nigra), stage 4 (temporal mesocortex and hippocampus), and stage 5 (insula, cingulate cortex, and temporal neocortex) belong to our PD-ICA network. Braak hypothesized that the medial temporal lobe served as a beachhead for further propagation to the remaining cortex, and interestingly three of the best supratentorial propagators identified here were the parahippocampal gyrus, anterior medial temporal lobe, and hippocampus (Table 3).

26Another prediction of the network-spread model is that connectome hubs should be especially vulnerable to disease spread. Hubs are defined as brain regions that are highly connected using graph theoretical metrics such as degree or betweenness. In theory, hubs, being highly connected, should have greater exposure to a toxic agent that is spread trans-neuronally. In a recent meta-analysis, found that most neurodegenerative diseases, including AD and fronto-temporal dementia, did indeed target hubs. However, one salient exception was PD. One possibility is that the PD studies included in the meta-analysis were too small in scale or methodologically incapable of detecting the true extent of damage. It is intriguing that our study identified numerous hub regions in the PD-ICA network, including the medial temporal lobe, putamen, insula, occipital cortex, anterior cingulate, superior frontal gyrus, and middle frontal gyrus. One notable hub region absent from our PD-ICA network is the posterior cingulate gyrus/precuneus, an area typically affected in AD. It would be interesting to see if this area eventually develops atrophy as PD progresses, and whether this is associated with cognitive impairment.

27The PD-ICA network demonstrated a correlation between atrophy and dopamine denervation measured with SPECT, in both PD patients and age-matched control subjects. Atrophy and dopamine denervation also both correlated with age in the control group. In sum, the control subjects demonstrated age-related loss of dopamine innervation and atrophy in the PD-ICA network. This finding is also novel. It is known that healthy aging is associated with a progressive loss of dopamine neurotransmission, and that this may account for motor slowing and executive cognitive impairment that occurs with age. Indeed, subtle Parkinsonian signs such as stooped posture, bradykinesia, and reduced facial expression in the healthy elderly were associated with SN neuron loss at post-mortem. Our results extend these findings by showing that neurodegeneration of the extended dopamine network also occurs in healthy aging, albeit without attaining the severity of outright PD, possibly via loss of neurotrophism, or perhaps, due to low-level mitochondrial dysfunction or synucleinopathy. PD might then target the PD-ICA network due to a dual hit effect of pathology superimposed upon normal aging.

28Data used in this article were primarily obtained from the PPMI database. In addition, rsfMRI and DW-MRI in healthy subjects were used to generate normal human connectomes for investigation of the network propagation model.

29The PPMI is described at www.ppmi-info.org. PPMI is a public–private partnership funded by the Michael J Fox Foundation for Parkinson's Research and funding partners listed at www.ppmi-info.org/fundingpartners. It is an observational, multicenter longitudinal study designed to identify PD biomarkers. Each participating PPMI site received approval from a local research ethics committee before study initiation and obtained written informed consent from all subjects participating in the study.

30For this study, we used the initial visit 3T high-resolution T1-weighted MRI scans acquired from September 2013 to January 2014. MRI data were acquired in 16 centers participating in the PPMI project, using scanners from three different manufacturers (GE medical systems, Siemens, and Philips medical systems). The acquisition parameters are detailed in the data set Website: http://www.ppmi-info.org/wp-content/uploads/2015/03/PPMI-MRI-Operations-Manual-V7-0-20JAN2015-FINAL.pdf.

31We also made use of the following data for each participant: age, sex, disease duration, score on the Movement Disorder Society—UPDRS III while off medications, score on the Montreal Cognitive Assessment battery and striatal binding ratio (SBR) from SPECT measurements with the tracer [123I]FP-CIT, a measure of dopamine neuron terminal density. In total data from 355 subjects (237 PD patients and 118 age-matched controls) were used. Six subjects, five PD patients and one age-matched control, were excluded from analysis due to failure in MRI processing. Clinical characteristics are shown in Table 4.

32We acquired rsfMRI in 51 healthy, right-handed volunteers (mean age: 23.6 ± 5.9, range: 18–47, 32 men, 63%). None of the subjects reported a history of drug abuse, neurological or psychiatric disorder. The experimental protocol was reviewed and approved by Research Ethics Board of MNI. All subjects gave informed consent. Scans were acquired using a Siemens MAGNETOM Trio 3T MRI system at the MNI. High-resolution, T1-weighted, three-dimensional volume acquisition for anatomic localization (1-mm3 voxel size) and resting-state echoplanar T2*-weighted images with blood oxygenation level-dependent (BOLD) contrast (3.5-mm isotropic voxels, TE 30 ms, TR 2 s, flip angle 90°) were acquired from all participants. Each resting-state scan was 5-min long (150 vol). Subjects were instructed to rest quietly with eyes open.

33To obtain white matter connectivity maps of normal brain, we used the Illinois Institute of Technology Human Brain Atlas v.3 (http://www.nitrc.org/projects/iit2/). This atlas contains structural (T1) and high angular resolution DW-MRI data, and probabilistic gray matter maps of the adult human brain in MNI space. The atlas was generated from MRI data from 72 human subjects (42 females (59%): mean age 26.6 ± 4.8 years, range 20–39 years; 30 males: mean age 31.9 ± 4.9 years, range 22–40 years).

34Local change in tissue density was calculated using DBM. DBM consists in spatially transforming each MRI non-linearly to a stereotaxic template, and using the local deformation as a measure of tissue expansion or atrophy. There are other methods to detect population differences in brain structure. VBM measures local gray matter density by transforming the brain to stereotaxic space, segmenting the tissue into gray and white matter and cerebrospinal fluid, and performing spatial smoothing on the gray matter maps so that local image intensity reflects gray matter density. We decided against VBM as it does not preserve the entirety of the MRI data, and there is some suggestion that it is less sensitive than DBM to subcortical atrophy. Another approach is to measure cortical thickness from the MRI, but this would also make it impossible to detect subcortical changes.

35For DBM, we registered every brain non-linearly to the MNI152-2009c template (available at http://www.bic.mni.mcgill.ca/ServicesAtlases/ICBM152NLin2009) and computed the deformation applied at each voxel using the procedure explained in. Pre-processing of MRI included denoising using optimized non-local means filtering, correction for intensity inhomogeneity, and linear intensity scaling using histogram matching to the MNI152-2009c template.

36The resulting images were registered to MNI space using the MNI152-2009c template, in two steps: (1) A hierarchical nine-parameter linear registration was computed between native MRI images and the template by maximizing the cross correlation of intensities as the similarity measure. The resulting transformation was applied to the MR image to resample it onto a 1-mm3 voxel grid and bring it into MNI space. (2) A hierarchical non-linear registration was performed on the linearly resampled scan to align it with the MNI152-2009c template. The resulting non-linear transformation field was inverted to generate a map of the deformations in template (MNI) space for each subject. Quality control was performed on each individual data set: the brain mask, and linear and non-linear registrations were visually inspected, and data sets with faulty registration were discarded.

37After the registration procedure, for each MRI and for each position in the brain x = x1, x2, x3, we obtain a displacement value in each direction to generate a field of vectors: U(x) = (u1(x), u2(x), u3(x)); that is, during the registration procedure position, x is displaced to the new position x + U(x) in the template space. This shows how much each voxel was moved from the MNI152-2009c template to match the subject's brain. To estimate local atrophy, an extra step is needed. Since a completely uniform displacement results in no change of volume, we are interested in the derivative of the displacement in each direction, that is, the determinant of the local Jacobian matrix of displacement.

38This Jacobian matrix is estimated using a first order approximation:where δ is the voxel dimension along the x2 axis. To calculate the relative ratio of the local volume change, we use the determinant of the Jacobian matrix |J| minus 1, that is, |J| − 1. By performing this calculation at each voxel, we obtain a map of local relative volumetric difference between each subject image and the MNI152-2009c template, reported in the MNI152-2009c template Talairach-like coordinate system.

39We created a composite anatomical atlas of gray matter regions in the cerebral hemispheres, cerebellum, and midbrain. For the supratentorial regions, we used the Hammers atlas (Copyright Imperial College of Science, Technology and Medicine, Alexander Hammers and University College London 2011. All rights Reserved) excluding the brainstem, cerebellum, and white matter. For the cerebellum, we used a public-domain high-resolution digital cerebellar atlas. These two atlases do not have adequate segmentation of three midbrain structures, the SN, subthalamic nucleus, and red nucleus. We therefore manually segmented these three structures using the Display software tool (McConnell Brain Imaging Centre) and three sources of anatomical information: the high-resolution MRI template (T1-weighted ICBM 2009c template, resolution = 0.5 mm3), the BigBrain, a 20-micron resolution digital brain atlas in MNI space, and the brainstem anatomical atlas of. The three structures were manually drawn on the high-resolution ICBM 2009c template. We then confirmed the segmentation of these three regions with a recently developed subcortical atlas based on ultra high-field MRI. The composite atlas thus created contains 112 cortical and subcortical structures (Figure 4—figure supplement 1, Figure 4—source data 1) and excludes all brainstem regions caudal to the SN.

40We used ICA to extract patterns of deformation in PD patients and age-matched controls. ICA is a statistical method to decompose multivariate data into statistically independent components. In this case, we used ICA to decompose the deformation maps into spatially distinct subcomponents. ICA was performed with MELODIC (http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/MELODIC), a toolbox that is part of the FSL analysis package. All the individual DBM images (patients and controls combined) were concatenated to create a 4-D image in which the first three dimensions are the individual 3-D deformation maps and the fourth dimension consists of the subjects. MELODIC applies an initial principal component analysis-based dimension estimation to find the optimal number of independent components and then uses this number in the decomposition procedure to identify the independent spatial components. Each of these components consists of a vector of normalized DBM values (one per subject) as well as a corresponding 3-D spatial map. The spatial maps were converted to z-statistic images via normalized mixture-model fitting, and thresholded at z = 3. The ICA algorithm in MELODIC is sensitive to sparsely distributed (super-Gaussian) data as typically seen in fMRI. We confirmed that the DBM data used here possessed this super-Gaussian property (kurtosis >4).

41To identify regions showing greater atrophy in PD, the average DBM values for each ICA component and each subject were entered in an unpaired t-test comparing PD subjects and age-matched controls (Bonferroni corrected for multiple comparisons). Also, the DBM values from each component were correlated with the age of each subject to identify patterns of deformation associated with aging. Again, Bonferroni corrections were applied. The DBM values from the component(s) that demonstrated a statistically significant group difference were compared to two disease-related clinical measures using linear regression: SBR and UPDRS III.

42In order to test the theory that PD targets normal brain networks, we analyzed rsfMRI data from 51 young healthy individuals. The rsfMRI data were preprocessed using the Neuroimaging Analysis Kit (NIAK), to perform slice timing correction, rigid body motion correction, and removal of slow temporal drift using a high-pass filter with 0.01 Hz cut-off. Physiological noise was accounted for by including white matter and cerebrospinal fluid signals as covariates. In the next step, the mean motion-corrected volume of each subject's fMRI data was first linearly (6 parameters: 3 translations, 3 rotations) registered to the native individual T1 image and then non-linearly registered to the MNI152 non-linear template (http://www.bic.mni.mcgill.ca/ServicesAtlases/HomePage). All data were resampled and smoothed with a 6-mm Gaussian kernel. All fMRI time series further underwent level one auto-regression (AR1) temporal de-noising.

43The mask of the SN described above was used to extract the average BOLD time series of each individual scan. The time series from right and left SN were averaged. Seed-based functional connectivity analysis was performed as implemented in fmristat (http://www.math.mcgill.ca/keith/fmristat/). The SN average BOLD signal was entered as a regressor in the design matrix and its correlation calculated with all voxels in the brain. A mixed effects model was applied to generate a t-map for the group. This was thresholded based on random field theory to achieve a p value of 0.05 after correction for multiple comparisons.

44We quantified the spatial similarity between PD-ICA atrophy maps and intrinsic networks in the healthy brain. We compared the PD-ICA map to intrinsic brain networks obtained from three different sources: the functional connectivity map obtained from the SN seed region, and ICNs from two published sources. We also compared the PD-ICA network map to a map derived from a meta-analysis of fMRI experiments where subjects tracked the value of rewarding stimuli.

45First, the functional connectivity map obtained from the SN seed-based analysis was compared to the 30 different ICA maps from the deformation analysis. Second, we compared the PD atrophy map (PD-ICA network) to the 70 ICNs identified by Smith et al. based on ICA analysis of resting-state data in healthy subjects. We further compared the PD-ICA map to 100 ICNs identified from the HCP rsfMRI data. Then, we compared the 30 ICA maps from the deformation analysis to a map identified in a meta-analysis of value-related fMRI studies, reasoning that dopamine networks implicated in PD would be similar to value networks in healthy brains (i.e., regions where BOLD signal tracks the value of experimental stimuli). Spatial cross-correlation (Pearson) was used as the measure of similarity. We chose |r| > 0.25 as indicative of similarity between two spatial maps as this has been argued to guard against false positives in a similar comparison between ICA-derived spatial maps.

46We generated the normal connectome of our 112 region brain atlas (52 paired bilateral regions, 8 midline regions), in two different ways, with rsfMRI and with DW-MRI. The goal here was to test the theory that geodesic (synaptic) proximity to the SN in healthy brain would predict the distribution of regional atrophy in PD.

47The modified Hammers-Cerebellum-Brainstem atlas was used to extract the average BOLD time-series, after correction for physiological noise, for all regions for each of the 51 rsfMRI acquisitions. All 112 regions in the atlas were used as separate masks and average time series were extracted for each. The time series used were the residuals of the linear model after correction for physiological noise and head motion. As a result, we obtained the correlation between all region pairs, which gives a 112 × 112 connectivity matrix for each acquisition. Then, a common group-based matrix called Rgroup was calculated. To do so, we followed steps explained in: each individual subject's (j) correlation matrix was converted to a standard normal metric using the Fisher transformation

48Then, all the Fisher transformed results were averaged.where N is the number of subjects. Finally, the group Z correlation matrix result was converted back to the correlation space using an inverse Fisher transform

49The modified Hammers-Suit-Brainstem atlas was also used to extract the average deformation values (Jacobian) for each region for each PD and age-matched control subject. A t-test was performed in each region to assess the difference between the PD and control groups, corrected for age. These t-values are a measure of the difference between the two groups and the sign of the t-value shows the direction of the effect, with negative values chosen to mean greater deformation (atrophy) in PD.

50The relationship between functional connectivity of each brain region with the selected seed region (SN) and the t-value of the deformation in each region (PD minus age-matched controls) was investigated using correlation analysis.

51Probabilistic anatomical connectivity values between each pair of atlas regions were estimated using a fully automated fiber tractography algorithm and the intravoxel fiber orientation distribution functions from the IIT Human Brain Atlas v.3. A maximum of 500-mm trace length and a curvature threshold of ±90° were imposed as tracking parameters. Based on the resulting voxel-region connectivity maps, the anatomical connection probability (ACP) between any pair of regions i and j (ACPij = ACPji) was calculated as the maximum voxel-region connectivity value between both regions. The ACP measure reflects the degree of evidence supporting the existence of each hypothetical white matter connection, independently of the density/strength of this connection, and is thus a measure of low susceptibility to gross fiber degeneration related to the aging processes. Then, effective anatomical distances between each region i and all the other ROI were estimated as the lengths of the shortest paths (in terms of ACP values) linking that region to all the other regions.

52Then, as for the functionally derived connectome above, we tested whether the PD minus control atrophy value in each atlas region depended on the geodesic distance to the SN in the DW-MRI connectome. We calculated the non-linear correlation between the regional T value of the difference in Jacobian between patients and controls and the effective anatomical distance to the SN.