Mapping track density changes in nigrostriatal and extranigral pathways in Parkinson's disease.

In Parkinson's disease (PD) the demonstration of neuropathological disturbances in nigrostriatal and extranigral brain pathways using magnetic resonance imaging remains a challenge. Here, we applied a novel diffusion-weighted imaging approach-track density imaging (TDI). Twenty-seven non-demented Parkinson's patients (mean disease duration: 5 years, mean score on the Hoehn & Yahr scale=1.5) were compared with 26 elderly controls matched for age, sex, and education level. Track density images were created by sampling each subject's spatially normalized fiber tracks in 1mm isotropic intervals and counting the fibers that passed through each voxel. Whole-brain voxel-based analysis was performed and significance was assessed with permutation testing. Statistically significant increases in track density were found in the Parkinson's patients, relative to controls. Clusters were distributed in disease-relevant areas including motor, cognitive, and limbic networks. From the lower medulla to the diencephalon and striatum, clusters encompassed the known location of the locus coeruleus and pedunculopontine nucleus in the pons, and from the substantia nigra up to medial aspects of the posterior putamen, bilaterally. The results identified in brainstem and nigrostriatal pathways show a large overlap with the known distribution of neuropathological changes in non-demented PD patients. Our results also support an early involvement of limbic and cognitive networks in Parkinson's disease.

Introduction

Parkinsonism is clinically defined by the presence of motor slowness with muscle rigidity, and/or tremor, and/or postural instability (Gibb, 1988). The most common cause of degenerative Parkinsonism in adults is Parkinson's disease. The main motor features of Parkinson's disease stem from the death of pigmented neurons in the substantia nigra (SN). Neuronal loss begins in the ventrolateral SN and this remains the most severely affected region throughout the illness (Damier et al., 1999a, Fearnley and Lees, 1991). As the disease progresses, pathology extends into the dorsal tier of the nigra and ventral tegmental area (Fearnley and Lees, 1994). Post-mortem studies have shown that the loss of cells in these regions results in profound dopamine (DA) depletion in the motor region of the striatum (Kish et al., 1988), with nigral projections to the dorsal and caudal putamen being most affected. A post-mortem diagnosis of Parkinson's disease requires evidence of cell loss in the SN, as well as the presence of Lewy pathology—a term that describes variously shaped insoluble intraneuronal deposits that contain misfolded α-synuclein. These aggregated α-synuclein deposits accumulate in neurites (Lewy neurites) and in neuronal somata (Lewy bodies). Furthermore, neuropathological abnormalities can be observed in other brainstem nuclei such as the pedunculopontine nucleus, locus coeruleus, and raphe nucleus (Gesi et al., 2000, Halliday et al., 1990, Jellinger, 1991, Zweig et al., 1993), and may extend up to the basal forebrain and cerebral cortex, possibly according to a predictable spatiotemporal pattern (Braak et al., 2003a, Braak et al., 2003b).

In vivo nigral and extranigral imaging biomarkers of Parkinson's disease are highly sought after for many reasons. Candidate biomarkers may aid premotor diagnosis and help differentiate Parkinson's disease from look-alike conditions such as essential tremor and atypical Parkinsonian syndromes. Perhaps more importantly, reliable imaging biomarkers may aid the development of disease-modifying therapies, as they can be used to monitor disease progression.

To date, no in vivo brain imaging modality has adequately captured the widespread spatial spectrum of brain abnormalities in non-demented patients with Parkinson's disease. Because a biochemical hallmark of Parkinson's disease is a deficiency of striatal DA, many imaging studies have focused on studying the problem directly. Measurement of striatal [18F]-l-dihydroxyphenylalanine (18F-DOPA) uptake with positron emission tomography (PET) is regarded by many as the “gold standard” for diagnosis of Parkinson's disease. However, while decreased 18F-DOPA uptake may also be observed in other brainstem regions and cortical areas (Pavese et al., 2012), it is generally acknowledged that it does not reveal the whole range of extranigral pathological abnormalities.

Despite recent advances in brain imaging, extranigral abnormalities remain difficult to capture in vivo in non-demented Parkinson's patients. As routine MRI is typically unremarkable in Parkinson's disease, the use of more advanced techniques is warranted. In this paper, we used track density imaging, an advanced diffusion-weighted imaging (DWI) method that allows the mapping of cerebral fiber pathways at a spatial resolution smaller than the voxel size of the original MRI (Calamante et al., 2010).

Diffusion-weighted imaging allows for the quantification of water mobility within tissue. In DWI analysis the movement of water within each voxel is modeled (e.g. using a tensor) and used for further analysis. To date, most studies have relied on scalar measures derived from the diffusion tensor model, such as fractional anisotropy (FA) and mean diffusivity (MD), which quantify the degree of anisotropy and average magnitude of local water diffusion, respectively. Most previous DWI studies in Parkinson's disease used FA and MD on focused regions of interest defined a priori, usually in the substantia nigra (Chan et al., 2007, Vaillancourt et al., 2009, Yoshikawa et al., 2004, Zhan et al., 2011). In general, results show decreased FA in Parkinson's patients compared to controls. Studies in rat (Soria et al., 2011) and mouse (Boska et al., 2007) models of Parkinson's disease have also reported similar decreases. Other studies using methods like tract-based spatial statistics, voxel-based FA analysis, and ROIs outside the substantia nigra, have resulted in less consistent reports. Changes in tensor-derived measures have been reported in the gyrus rectus (Ibarretxe-Bilbao et al., 2010), the genu of the corpus callosum, the superior longitudinal fasciculus (Gattellaro et al., 2009), as well as motor and frontal cortices (Zhan et al., 2011). As these past studies have relied on DWI sequences with relatively low numbers of unique gradient directions [6 to 64, (Cochrane and Ebmeier, 2013)], they may have lacked the sensitivity to show disturbances using unbiased voxel-based analyses.

Here we acquired data in 27 non-demented Parkinson's disease patients and 26 matched controls using an advanced DWI sequence, in which 120 diffusion gradients were applied, allowing for the fitting of higher order models of the signal in each voxel. We reconstructed white matter streamlines from DWI images with constrained spherical deconvolution, which generally provides far greater accuracy than alternatives like single or multi-tensor models (Tournier et al., 2004, Tournier et al., 2008). In this paper, when referring to the results of fiber tractography, we will use the terms “fibers,” “tracks,” and “streamlines” interchangeably, and these should not be confused with actual biological tracts (Jones et al., 2013).

The added value of track-weighted imaging comes from the information obtained by tracking the neural pathways (Calamante et al., 2012b). Track density imaging in essence is just the resampling of the fiber track data into user-specified volumetric data (Calamante et al., 2010). Sampling the tractography dataset with a voxel size smaller than the original diffusion-weighted scan is known as “super-resolution” track-weighted imaging. In the approach presented here – track density imaging (TDI) – the manufactured signal is simply the number of streamlines passing through each voxel (Calamante et al., 2010). A similar approach, average pathlength mapping, has shown promise for evaluating traumatic brain injuries (Pannek et al., 2011). Streamline density imaging has been validated histologically for its ability to resolve white matter boundaries (Calamante et al., 2012a).

The present study also differed from most previous DWI studies by using a whole-brain voxel-based approach, rather than regions-of-interest for data analysis. This allowed us to explore and test for changes in both nigral and extranigral areas in the patients as compared with the control groups. Based on the spatial distribution of neuropathological abnormalities reported in non-demented Parkinson disease patients, we hypothesized that differences in track density would appear predominantly in the brainstem and nigrostriatal pathways.

Materials and methods

Participants

We studied 27 (14 males) patients clinically diagnosed with Parkinson's disease and 26 (14 males) healthy control subjects from a larger data sample after excluding those with poor quality imaging data. The groups were matched by age, sex, and highest achieved education level. Disease stage within the patient population was assessed in the “on” state using the Hoehn & Yahr scale (Hoehn and Yahr, 1967). Disease severity was evaluated using the Unified Parkinson's Disease Rating Scale (UPDRS) (Fahn and Elton, 1987). Quality of life was estimated using the Parkinson's Disease Questionnaire (PDQ39) (Jenkinson et al., 1997). Subjects were also administered several psychological tests: global cognitive function [Mattis Dementia Rating Scale (Mattis, 1988), MMSE (Folstein et al., 1975), Symbol Digit Modalities test, verbal fluency test], inhibition [Stroop test, random number generation (Jahanshahi et al., 2006)], episodic memory [Rey auditory verbal learning test (Rey, 1958)], updating of working memory [letter running span memory task], cognitive flexibility [Modified Wisconsin Card Sorting Test], visuospatial judgment [Judgment of Line Orientation (Benton et al., 1983)], and anxiety [Hospital Anxiety and Depression Scale (Zigmond and Snaith, 1983)].

Twenty-four of the 27 patients were taking a combination of several classes of drugs: levodopa (immediate and controlled release), nonergot-derived dopamine receptor agonists (Pramipexole, Ropinirole), and a monoamine oxidase B inhibitor (Rasagiline). The remaining 3 patients were not taking any anti-Parkinsonian medications at the time of scanning. Levodopa and dopamine agonist dosages were pooled and summarized as the levodopa equivalent daily dose [LEDD, (Hobson et al., 2002)]. Total daily l-DOPA equivalent dosages ranged from 0 to 900 mg. Demographic and clinical data are summarized in Table 1. Written informed consent was obtained from all participants in accordance with the Declaration of Helsinki. The Ethics Committee of the University of Liège approved the study.

Table 1

Demographics for the study cohort.

	Healthy controls (n = 26)	Parkinson's patients (n = 27)	t-Test, p-value
Age	64 (8)	66 (8)	0.549
Sex (M:F)	14:12	14:13
Years of education	13 (3)	11 (3)	0.133
ICV (mm3)	1478 (152)	1516 (148)	0.360
BMI (kg/m2)	25 (3)	25 (3)	0.660
Hand dominance (L:R)	2:24	2:25
Hoehn & Yahr stage		1.5 (0.62)
Most affected side (L:R)		10:17
Disease duration (years)		5 (3)
LEDD (mg)		323 (255)
UPDRS Section 2		9 (6)
UPDRS Section 3		14 (7)
Mattis	139 (4)	136 (4)	0.004
MMSE	29 (1)	28 (1)	0.022
HADS total	10 (4)	13 (6)	0.066
PDQ39 mobility		20 (18)
PDQ39 total		189 (114)
Rey Auditory Verbal Learning Test	53 (11)	44 (11)	0.006
SDMT	51 (10)	45 (12)	0.062
JOLO	27 (4)	25 (4)	0.055

Values reflect mean (st. dev.). Two-tail t-tests were performed with an assumption of unequal variance in each group. ICV = intracranial volume, BMI = body mass index, LEDD = l-DOPA equivalent daily dose (Hobson et al., 2002), UPDRS = Unified Parkinson's Disease Rating Scale, MMSE = Mini Mental State Examination, HADS = Hospital Anxiety and Depression Scale, PDQ = Parkinson's Disease Questionnaire, SDMT = Symbol Digit Modalities Test, JOLO = Judgment of line orientation test.

Imaging data acquisition

Data were acquired on a 3 T head-only MR scanner (Magnetom Allegra, Siemens Medical Solutions, Erlangen, Germany) operated with an 8-channel head coil. Diffusion-weighted (DW) images were acquired with a twice-refocused spin-echo sequence with EPI readout at two distinct b-values (b = 1000, b = 2500 s/mm2) along 120 encoding gradients that were uniformly distributed in space by an electrostatic repulsion approach (Jones et al., 1999). This sequence is designed specifically to reduce the distortions induced by eddy-currents in the diffusion-weighted images (Reese et al., 2002). For the purposes of motion correction, 22 unweighted (b = 0) volumes, interleaved with the DW images, were acquired. Volumes were acquired with a repetition time (TR) of 6800 ms, an echo time (TE) of 91 ms, and a field-of-view (FOV) of 211 mm2. Maximum slew rate was 400 mT/(m/ms) and maximum gradient amplitude was 40 mT/m. No parallel imaging techniques were used. The multi-channel head coil was used to increase the signal-to-noise ratio, and not the speed of the acquisition. Volumes were acquired with a 6/8 partial Fourier factor. Voxels were isotropic with dimensions of 2.4 × 2.4 × 2.4 mm3 and volumes were acquired in 54 transverse slices using an 88 × 88 voxel matrix. Total acquisition time was approximately 35 min. MRI data were acquired while patients were on their usual anti-Parkinsonian medication so they would remain as still as possible in the scanner.

Structural and quantitative maps of T1, T2*, proton density (PD) and magnetization transfer (MT) at 1 × 1 × 1 mm3 resolution were calculated from a multi-parameter protocol based on a 3D multi-echo fast low angle shot (FLASH) sequence (Weiskopf and Helms, 2008). Three co-localized 3D multi-echo FLASH datasets were acquired with predominantly proton density weighting (PDw: TR/α = 23.7 ms/6°), T1 weighting (T1w: 18.7 ms/20°), and MT weighting (MTw: 23.7 ms/6°; excitation preceded by an off-resonance Gaussian MT pulse of 4 ms duration, 220° nominal flip angle, 2 kHz frequency offset). The PDw acquisition used eight bipolar gradient echoes at equidistant TE between 2.2 ms and 19.7 ms, whereas the T1w and MTw acquisition used six bipolar gradient echoes at equidistant TE between 2.2 ms and 14.7 ms. All three had 425 Hz/pixel bandwidth to avoid susceptibility-related distortions. Volumes were acquired in 176 sagittal slices using a 256 × 224 voxel matrix. The total acquisition time was approximately 19 min.

A B1 map was calculated using the actual-flip-angle imaging method (Yarnykh, 2007) based on two interleaved 3D FLASH acquisitions (repetition times TR1 = 33 ms, TR2 = 165 ms, TE1 = TE2 = 3.05 ms, nominal flip angle = 60°, acquisition time = 5 min.). The B1 map was used to correct the multiparameter maps for B1 field bias (Volz et al., 2010).

Visual inspection of the raw images was performed to ensure that subject data had no (i) ghosting, (ii) signal dropout, (iii) major motion distortions, or (iv) scanner-induced vibration artifacts (Tournier et al., 2011).

Imaging data processing

Interleaved unweighted images from the diffusion sequence were realigned to the first unweighted volume with a rigid body transformation using SPM8 (Wellcome Trust Centre for Neuroimaging, UCL, UK). Registration was performed (rigid, mutual information) between the first unweighted volume and each of the interleaved unweighted volumes (e.g. register b03 to b01). The translation and rotation values between b01 and b022 were linearly interpolated and applied to the weighted volumes. The assumption behind this procedure is that the subject's displacement was linearly continuous in between unweighted volumes. This put all the weighted images in alignment with the first b0 volume without the contrast problems of co-registering weighted and unweighted images. Diffusion gradient vectors were rotated accordingly (Leemans and Jones, 2009). For each diffusion-weighted volume, a non-local mean filter was applied (Maggioni et al., 2013) and noise was corrected using power image correction adapted for multi-coil acquisitions (André et al., 2014). No further corrections were applied to correct for eddy current-induced distortions in the diffusion-weighted volumes because the diffusion sequence did a sufficient job of suppressing them.

We included motion as a nuisance regressor in our analysis by computing a version of the total motion index (Yendiki et al., 2014) from the rigid transformations between unweighted volumes. We calculated average volume-by-volume translation between unweighted images by averaging the magnitude of the translation vectors between each successive unweighted image. Similarly, average volume-by-volume rotation was computed by averaging the sum of the absolute values of the rotation values between successive unweighted images.

FA and MD maps

A brain mask was generated using the “median Otsu” method as implemented in Dipy (version 0.7.1) (Garyfallidis et al., 2014). In this method, a median filter is applied to smooth the unweighted volume, and then an automated histogram approach is used to separate the foreground and background. Using the low b-value (b = 1000 s/mm2) data, tensors were fit at each voxel using non-linear least squares, and fractional anisotropy and mean diffusivity maps were also generated from the tensors.

TDI maps

The fiber response model was estimated for each subject from the high b-value (b = 2500 s/mm2) diffusion-weighted images. A mask of single fiber voxels was extracted from the thresholded and eroded FA images. Only strongly anisotropic (FA > 0.7) voxels were used to estimate the spherical-harmonic coefficients of the response function (Tournier et al., 2004, Tournier et al., 2008). Using non-negativity constrained spherical deconvolution, fiber orientation distribution (FOD) functions were obtained at each voxel.

For both the response estimation and spherical deconvolution steps we chose a maximum harmonic order of 8. Seed masks were created from white matter probability maps, estimated from signal variation in the DWIs,2 thresholded above 0.4. This threshold was chosen by trial and error, and is in general lower than what would be used on young healthy subjects, possibly due to variations in levels of iron (Graham et al., 2000). Probabilistic tractography was performed using randomly placed seeds within subject-specific white matter masks. Fiber tracking settings were as follows: number of tracks = 5,000,000, FOD magnitude cutoff for terminating tracks = 0.1, minimum track length = 10 mm, maximum track length = 200 mm, minimum radius of curvature = 1 mm, tracking algorithm step size = 0.2 mm. Streamlines were terminated when they (i) extended out of the white matter mask, or (ii) could not progress along a direction with an FOD magnitude or curvature radius higher than the minimum cutoffs.

Spatial normalization of the fiber tracks followed a procedure laid out previously (Pannek et al., 2011). First, track density and total pathlength maps (TPM) were generated with 1 × 1 × 1 mm3 isotropic voxels, and used to create average pathlength maps (TPM/TDI = APM) for each subject. The APMs were then rigidly aligned to the MNI152 1 mm T1 template with mutual information and averaged to create an initial template. The APMs were then used to build a template through 4 iterations, beginning with an initial affine transformation followed by greedy symmetric diffeomorphic normalization (SyN) with cross-correlation (Avants et al., 2010, Avants et al., 2011). For each subject, the transformations (rigid, affine, warp) were then individually inverted and applied in inverse order to a unit warp field generated in the final template space.3 This was then used to normalize the tracks into template space with MRtrix, by warping each point along the fibers. Finally, track density images were generated from the normalized tracks, again with 1 × 1 × 1 mm3 isotropic voxels. Normalization steps were performed using Advanced Normalization Tools (ANTS, http://www.picsl.upenn.edu/ANTS/). No spatial smoothing was applied.

Motion-corrected DWIs were assessed to ensure that the correction procedure performed adequately. Inspection of the scalar images, such as FA, MD, seed masks, single-fiber response estimation masks, APM, and pre- and post-normalization TDI was also performed. Track datasets were truncated to 100,000 tracks to assess registration with the template using MRview (Tournier et al., 2012). Diffusion-weighted image and fiber track processing was performed with the MRtrix package version 0.2.12 (J-D Tournier, Brain Research Institute, Melbourne, Australia, http://www.brain.org.au/software/) (Tournier et al., 2012). The tractography and spatial normalization workflows were developed in Python with the Nipype pipeline architecture and are freely available online (Gorgolewski et al., 2011).

MT maps

For display purposes, semi-quantitative MT saturation maps were calculated within SPM8. These MT maps show the percentage loss of magnetization imposed by a single MT pulse (Draganski et al., 2011). The contrast in these MT images, which is higher than in typical T1-weighted images, enables more accurate distinction between gray and white matter, especially for the basal ganglia and substantia nigra (Helms et al., 2009). MT maps were processed using unified segmentation in order to create masks of gray matter, white matter, and cerebrospinal fluid (Ashburner and Friston, 2005). Diffeomorphic registration (DARTEL) was used to normalize the MT maps to a common template (Ashburner, 2007). Target images for the DARTEL registration were modified tissue probability maps from previous MT segmentations (Draganski et al., 2011), which included detailed segmentation of subcortical regions.

To display the final results, the study mean MT map was coregistered to the average TDI map with an initial rigid and secondary affine transformation using mutual information. Greedy symmetric diffeomorphic normalization with Mattes mutual information was used register the mean TDI map with the MNI152 1 mm T1 template. The identified transformation was applied to the mean MT image. Registration was performed in ANTs and accuracy was assessed visually at each step.

Statistical analysis

FA and MD data

Voxelwise statistical analysis of the FA data was carried out using tract-based spatial statistics (TBSS) (Smith et al., 2006), part of FSL 5.1 (Smith et al., 2004). TBSS projects all subjects' FA data onto a mean FA tract skeleton, before applying voxelwise cross-subject statistics. TBSS was also carried out on the mean diffusivity (MD) images. Non-parametric statistics were performed using FSL's Randomise (5000 iterations), and 2D Threshold-Free Cluster Enhancement (TFCE) (Smith and Nichols, 2009) with the mean-centered total motion index as a nuisance covariate.

TDI data

An optimal white matter mask was generated from all of the track density images using the SPM8 Masking toolbox and the Luo-Nichols anti-mode method of automatic thresholding (Luo and Nichols, 2003, Ridgway et al., 2009). Permutation testing was performed using Randomise (5000 iterations) with TFCE and the mean-centered total motion index as a nuisance covariate. Only clusters with p < 0.05 after correcting for family-wise error rate were considered statistically significant. Tentative anatomical labels were estimated from a number of freely available brain atlases. Brainstem regions were examined by comparison with a high-resolution single subject DTI scan (Aggarwal et al., 2013). Cortical clusters were localized with the JHU ICBM-DTI-81 white-matter atlas (Oishi et al., 2010), the JHU white-matter tractography atlas (Hua et al., 2008), and the Jülich Histological Atlas (Bürgel et al., 2006).

Pearson correlation was used to explore interactions with numeric clinical metrics (e.g. UPDRS III score). The Spearman rank-order correlation coefficient was used for non-numeric variables (e.g. Hoehn and Yahr stage). Bonferroni correction was used to account for the number of independent statistical tests. Between-group differences in psychological metrics, as well as total motion index, were assessed with a two-sample t-test with an assumption of unequal variance.

Data were processed using the NITRC Computational Environment for Cluster Compute Instances using Amazon Elastic Compute Cloud (EC2) servers. The virtual machine snapshot is available upon request.

Results

Behavioral data

In our cohort the Parkinson's disease patients scored significantly lower than controls on the Mattis scale, mini-mental state examination, and the episodic memory test (Rey auditory verbal learning test), suggesting mild cognitive disturbances. Apart from these, there were no statistically significant differences between the patients and controls on any of the psychological measures. Every patient's performance was above the standard cutoff threshold for the Mattis dementia rating scale (Llebaria et al., 2008), and none of them met standard criteria for dementia associated with Parkinson's disease (Emre et al., 2007). The group mean total motion index was larger in PD patients than controls, but this did not reach statistical significance (p = 0.053).

Tract-based spatial statistics (TBSS) on FA and MD maps

No statistically significant results were obtained from the TBSS analysis for either FA or MD images.

Voxelwise analysis of TDI map

Statistically significant increases in track density in Parkinson's disease, relative to controls, were found in brainstem and extrapyramidal motor networks, limbic, and cognitive circuits. Areas of increased TDI were strikingly symmetric between hemispheres. Fig. 1 shows axial slices of the population mean MT image, with the clusters containing significantly increased TDI overlaid. There were no statistically significant decreases in track density in the Parkinson's patients, compared to controls. To clarify the widely distributed but significant results, we lowered the statistical threshold to 0.01 (family-wise error rate corrected). The percent signal change in PD for clusters significant at pFWE < 0.01 is documented in Table 2.

Fig. 1

Overview of the increases in track density. Group-level statistical map of increased streamline density throughout the brain. Clusters are overlaid on the group-mean magnetization transfer image. MNI coordinates in millimeters are displayed on top of each slice. Clusters shown in blue are significant at pFWE < 0.05 and those in red are significant below pFWE < 0.01, estimated with threshold-free cluster enhancement.

Table 2

Cluster details.

Size (mm3)	Mean number of tracks per voxel			Location	Side	Centroid (mm)	Peak (mm)
	Controls	Parkinson's disease	Percent signal change in PD
3062	416.5	708.6	70%	Brainstem, substantia nigra, superior cerebellar peduncle, nigrostriatal pathway	Bilateral	0, − 31, − 26	− 2, − 40, − 49
1267	511.3	733.0	43%	Optic radiation	Right	32, − 72, 16	39, − 70, 18
1030	674.5	981.7	46%	Posterior inferior fronto-occipital fasciculus	Right	41, − 48, 9	38, − 45, 11
814	685.6	959.6	40%	Superior longitudinal fasciculus, parietal lobe	Right	36, − 18, 31	37, − 24, 29
434	324.2	517.2	60%	Middle cerebellar peduncle	Bilateral	− 3, − 23, − 39	6, − 17, − 39
302	477.2	697.4	46%	Anterior inferior longitudinal fasciculus	Right	44, − 23, − 15	46, − 14, − 25
133	114.2	246.4	116%	Lateral occipital cortex	Right	55, − 62, 1	59, − 62, − 1
109	935.9	1474.1	58%	Superior longitudinal fasciculus, parietal lobe	Left	− 33, − 20, 27	− 32, − 18, 26
97	539.6	767.6	42%	Superior longitudinal fasciculus, parietal lobe	Right	42, − 36, 34	45, − 37, 35
85	816.3	1218.7	49%	Superior longitudinal fasciculus, parietal lobe	Left	− 29, − 23, 37	− 30, − 21, 37
75	1112.5	1684.6	51%	Superior longitudinal fasciculus, parietal lobe	Right	35, − 39, 23	36, − 40, 23
38	231.9	433.8	87%	Posterior cingulum	Left	− 11, − 36, 37	− 11, − 36, 36
34	713.7	1035.9	45%	Corona radiata near somatosensory cortex	Left	− 20, − 44, 43	− 20, − 44, 44
25	370.6	544.7	47%	Anterior superior longitudinal fasciculus	Left	− 41, − 13, 28	− 41, − 12, 27
20	657.9	1044.8	59%	External capsule	Left	− 32, − 4, − 12	− 33, − 4, − 12
18	538.6	790.9	47%	External capsule	Right	30, − 13, − 9	30, − 14, − 9
17	982.3	1478.5	51%	Posterior optic radiation	Right	32, − 56, 18	32, − 55, 18
14	206.1	320.5	56%	Occipital cortex	Right	29, − 79, 27	29, − 78, 26
10	318.5	604.1	90%	Precuneus	Right	16, − 42, 43	16, − 42, 43

Information about significant (pFWE < 0.01) clusters identified in the study. Percent signal change in Parkinson's disease within each cluster was calculated using (TDIPD − TDIHC)/TDIHC. Clusters were largely bilateral. The largest cluster encompassed the substantia nigra and extended upwards to the striatum, bilaterally.

Brainstem and extrapyramidal motor networks

The brainstem and cerebellar peduncles showed widely increased TDI in the patients. Results projected onto the brainstem portion of the superior, middle, and inferior cerebellar peduncles, as well as the pontine crossing tract and the medial lemnisci. Brainstem nuclei also affected were the locus coeruleus, pedunculopontine nucleus, ventral tegmental area, and raphe nucleus.

Increased track density was found bilaterally in a pattern appearing to connect the substantia nigra and striatum. The significant regions encompassed the caudal ventral substantia nigra and extended dorsally toward the medial aspects of the posterior putamen. A coronal slice through the substantia nigra with the significantly (pFWE < 0.01) affected areas overlaid is shown in Fig. 2.

Fig. 2

Substantia nigra cluster. Single coronal slice showing the cluster that encompasses the substantia nigra (pFWE < 0.01). The cluster is overlaid on the study mean MT map, so that the SN is more visible. This cluster extends upward through the cerebral peduncles and terminates at the posterior medial putamen.

Limbic circuits

TDI increases were found in the ventral tegmental area and in white matter regions surrounding the ventral striatum both medially and laterally. The pallidum, as well as the anterior and dorsomedial nuclei of the thalamus were also affected. At the cortical level, the entire cingulum was affected heavily in the right hemisphere, but only in posterior regions of the left hemisphere. Bilateral increases in TDI were also found in the orbitofrontal cortex. Fig. 3 shows slices of affected areas in the limbic system.

Fig. 3

Limbic system clusters. Limbic involvement was widespread, involving the cingulum, orbitofrontal cortex, anterior and dorsomedial thalamus. Clusters are overlaid on the study mean TDI map. Clusters in blue are significant below pFWE < 0.05 and clusters in red are significant at pFWE < 0.01.

Cognitive circuits

Increased streamline density in Parkinson's patients was found in many areas of the cerebral white matter. Major associative pathways such as the superior and inferior longitudinal fasciculi (SLF, ILF) showed sections of increased TDI.

In the frontal cortex, the genu of the corpus callosum and areas of the anterior corona radiata were affected. The superior parietal lobule was widely affected, ostensibly by changes to the superior longitudinal fasciculi. Regions in the primary somatosensory cortices (BA2, BA3b) were notably affected. The anterior temporal poles were clearly affected by variations along the inferior longitudinal fasciculi. The occipital cortex was also affected, specifically in regions of the inferior fronto-occipital fasciculi, posterior thalamic radiations, optic radiations, and the forceps major. Fig. 4 shows affected areas involved in cognition.

Fig. 4

Cognitive network clusters. White matter tracts involved in cognitive function were notably affected. The superior and inferior longitudinal fasciculi, superior parietal lobule, and occipital cortex showed increased track density. Clusters are overlaid on the study mean TDI map. Clusters shown are significant at pFWE < 0.01.

Discussion

This whole brain voxel-based advanced MRI study identified brainstem and nigrostriatal pathways as significantly affected brain areas in non-demented Parkinson's disease patients as compared with controls. To the best of our knowledge, this is the first time a single MRI method has identified abnormalities projecting onto these areas. This is consistent with the known distribution of neuropathological abnormalities in early-stage Parkinson's disease (Braak et al., 2003b).

Our whole-brain voxel-based analysis also revealed significantly increased streamline density in cognitive and limbic networks. Since patients enrolled in the present study were in the first years after the diagnosis, dementia-free, and without major cognitive or behavioral deficits at the time of scanning (Table 1), it would be valuable to test if these imaging findings are predictive of further clinical deterioration.

Brainstem and extrapyramidal motor circuits

Dopaminergic denervation within the nigrostriatal pathway is the most well described neuropathological feature of Parkinson's disease, but to date, no MRI modality has been sensitive enough to clearly identify alterations to this pathway. The method used here was able to capture microstructural abnormalities within both the SN and fiber tracts to the striatum. Track density abnormalities predominated in the ventrolateral portion of the substantia nigra. This is consistent with neuropathological observations in Parkinson's disease (Fearnley and Lees, 1991), as the greatest dopaminergic cell loss occurs in nigrosome 1 of the substantia nigra pars compacta, which is in the caudal and ventrolateral region of the SN (Damier et al., 1999b).

Compared to controls, patients showed increased fiber density from the SN to the posterior part of the putamen, which is predominantly involved in motor circuits with the (pre)motor cortices. This finding is in agreement with post-mortem studies in Parkinson's patients that typically show profound dopamine depletion in the motor region of the striatum (Kish et al., 1988), with nigral terminals in the dorsal and caudal putamen being most affected.

We were unable to replicate findings of decreased FA at the level of the substantia nigra (Chan et al., 2007, Peran et al., 2010, Vaillancourt et al., 2009, Yoshikawa et al., 2004, Zhan et al., 2011). There are many possible factors that could explain the lack of significant changes in tensor-derived indices in our sample. First, DWI scanning parameters differ enormously between studies, with many of the previously published studies reporting data with 1.5 T scanners, using thick slices (which can include peri-nigral fiber tracts), and applying very few diffusion directions. Second, nearly all of these studies use ROIs drawn on FA or T1-weighted images, neither of which clearly delineates the substantia nigra. Third, the number of subjects is also quite low in many studies, and patient disease stage varies both within and between studies. Furthermore, our patient population is fairly early stage (mean Hoehn and Yahr stage: 1.5), and though pathological damage to the SN occurs well before symptom onset, it also continues as the disease progresses. A recent meta-analysis of the literature concluded that there is no significant disease effect on mean diffusivity (Schwarz et al., 2013), and that there is extremely large variation in the results presented on FA. They also reported that after excluding studies with unusually high values of nigral FA in the control groups, there was no disease effect on FA, either.

Increased TDI was found in caudal brainstem areas mainly in its posterior aspect. We must acknowledge that the spatial resolution of diffusion weighted MRI might be insufficient for accurate labeling of results given the anatomical complexity of the brainstem. Anatomical labeling of the involved structures is tentative, rather than definitive, and should be interpreted with caution. Significant clusters projected onto the known location of the locus coeruleus and pedunculopontine nucleus, ventral tegmental area, raphe nucleus, and pontine crossing tract.

The pedunculopontine area is heterogeneous in composition and contains cholinergic, GABAergic, and glutamatergic neurons. The PPN is heavily affected in Parkinson's disease, and the loss of roughly 50% of its cholinergic neurons (Jellinger, 1988) has been correlated with disease stage (Rinne et al., 2008). Neuropathological abnormalities in the pedunculopontine area and its connections with basal ganglia, cerebral cortex, cerebellum, and spinal cord are considered to be pivotal, not only for disturbances of gait and posture, but also for non-motor problems related to arousal and cognition (Benarroch, 2013, Stefani et al., 2013).

These results are also consistent with the reduction in medium-sized locus coeruleus neurons found in Parkinson's disease (Patt and Gerhard, 1993). In the post-mortem brain tissue of Parkinson's patients the average loss of noradrenergic neurons compared to controls in this region is around 70% (Bertrand et al., 1997, Zarow et al., 2003). Furthermore, connections between the locus coeruleus and the cerebellum, thalamus, and cerebral cortex are known to influence locomotion (Watson and McElligott, 1984) and arousal. Detriments in both of these systems are considered cardinal signs of Parkinson's disease (Lima, 2013).

In this study the proximal parts of middle and superior cerebellar peduncles were affected bilaterally. Previous ROI-based studies did not identify any changes in these areas in Parkinson's disease (Nicoletti et al., 2008, Paviour et al., 2007), although both relied on relatively poor quality diffusion data. The superior cerebellar peduncles contain fibers that connect the locus coeruleus and pedunculopontine nucleus with deep cerebellar nuclei and the cortex (Hazrati and Parent, 1992, Liu, 2012). Clinically relevant cerebellar involvement is uncommon in Parkinson's disease, though, and methodological issues may contribute collectively to the results identified in this region.

In the raphe nucleus, evidence for neuronal loss is less consistent. Some post-mortem studies of patients with Parkinson's disease have revealed a profound loss of 5-HT neurons (Jellinger, 1991) and their terminals (Birkmayer and Birkmayer, 1987, Kish et al., 2008), although other studies have not. This variability may contribute to the heterogeneity in mood disturbances found in Parkinson's patients.

Limbic network

Increased TDI was found bilaterally in limbic areas including the ventral tegmental area, ventral striatum, medial thalamus, orbitofrontal and anterior cingulate cortices.

The VTA loses considerable amounts of dopaminergic neurons in patients with Parkinson's disease (Uhl et al., 1985). The ventral tegmentum is populated largely by dopaminergic neurons and is tightly linked with the reward circuit of the nucleus accumbens (Sesack and Grace, 2009). It projects to the forebrain and limbic system and a loss of integrity in these pathways may influence patients' emotional state. In this study the patients showed a trend toward increased anxiety and depression, but scores on the HADS scale did not reach the threshold for the diagnosis of depression (Table 1).

Regions corresponding to the anterior and dorsomedial thalamic nuclei were found to have increased TDI in PD. The anterior, dorsomedial, and ventral anterior thalamic nuclei have previously been shown to have reduced FA in Parkinson's disease (Planetta et al., 2013). Both the anterior and dorsomedial nuclei are involved in the limbic system. The anterior nucleus projects to the cingulate gyrus, whereas the dorsomedial nucleus connects to the amygdala, anterior cingulate, and prefrontal cortices.

The patients in this study were in early stages of the disease and showed no evidence of dementia. For this reason we did not expect to find substantial abnormalities in the cerebral cortical areas. Neuropsychological testing, however, did indicate very mild cognitive disturbances in the patients.

Clusters of significantly increased TDI in patients projected primarily onto posterior areas of the cerebral white matter. The corresponding cortical areas consistently show disturbed resting-state activity in PD populations studied using functional imaging methods such as 18F-fluorodeoxyglucose PET and arterial spin labeling (Garraux et al., 2011, Hu et al., 2000, Melzer et al., 2011, Tang et al., 2010). Our study supports the hypothesis that altered activity in these areas may be related to damage within white matter pathways targeting these areas. This may contribute to the mild cognitive disturbances found in PD patients when compared with controls. Our transversal study design does not allow us to test whether or not altered white matter tracts in posterior areas are predictive of future cognitive deterioration.

Methodological considerations

Our analysis, unexpectedly, showed an increase in streamline density in the PD patients. The reason for this increase is unclear, though, so we are forced to speculate as to its cause. We have identified several potential sources of error. The first relates to the calibration and use of constrained spherical deconvolution to fit fiber orientation distribution models to the DWI data.

The key assumption behind spherical deconvolution approaches is that the measured signal in any voxel is a convolution of the fiber orientation distribution in that voxel and the diffusion-weighted signal response for a single-fiber population. The single-fiber response in this study is estimated on a subject-by-subject basis and assumed to be constant throughout the brain. The fiber response functions are approximated by a linear combination of spherical harmonic basis functions, which in this study, were truncated at the 8th order harmonic. The impact of higher harmonic orders (than 8) on the spherical harmonic fit has recently been shown to be negligible for the b-value used in this study (b = 2500 s/mm2) (Tournier et al., 2013). In this study, highly anisotropic (FA > 0.7) voxels were assumed to contain a fiber population oriented in a single direction, and used for fiber response calibration. These were extracted from an eroded FA image so that voxels at the edges of the brain, which are more prone to artifacts, were not considered. Future studies may consider adding a calculation of tensor mode (Ennis and Kindlmann, 2006) to help select voxels for calibration. Iterative (Tournier et al., 2013) and recursive (Tax et al., 2014) methods for CSD calibration have also recently been demonstrated. In our study the calibration values were not found to differ between groups (two sample t-test, unequal variance, data not shown).

One possible source of error is potential incongruity between the approximated fiber response function and the actual fiber response function found in chronically damaged fiber tracts. It has been demonstrated that poor estimates of fiber response function, (e.g. caused by calibrating the response function from voxels with multiple fiber populations) can lead to spurious fiber orientations being resolved (Parker et al., 2013). These artifactual fiber orientations may be especially apparent in regions with complex micro- and macro-structural properties, multiple fiber populations, or partial volume effects. We applied the standard threshold of FOD magnitude > 0.1, which substantially reduces the number of spurious peaks capable of disrupting the tractography algorithm (Parker et al., 2013).

It remains unclear how researchers should proceed when performing spherical deconvolution within pathological populations containing local changes in microstructural properties. Image-wide fiber response calibration may not be the best solution for pathological brains. Lowering the maximum harmonic order allowed in constrained spherical harmonic deconvolution (CSD) can reduce its sensitivity to miscalibration (Parker et al., 2013). Other alternatives to CSD that may reduce spurious fiber orientations are damped Richardson–Lucy deconvolution (Dell'Acqua et al., 2010) and voxel-by-voxel response estimation (Anderson, 2005).

The second concern involves the white matter masking and tractography seeding used in this study. We estimated white matter masks from the signal variation in the DWIs, and because this technique was used, they may have differed in spatial extent between subjects. White matter masks should be segmented carefully from high-resolution structural images and subsequently co-registered to diffusion space. It may also be beneficial to seed a specific number of tracks from every white matter voxel, rather than a large number from randomly chosen masked voxels.

One notable subjective choice in this study was the resolution of the TDI and APM maps. We chose 1 mm3 isotropic voxels based on our own testing and because track quantification studies at this resolution have been previously published (Besseling et al., 2012, Pannek et al., 2011). The primary requirement for the TDI and APM resolution is that it should provide reliable maps with small variability across subjects, so that template building and normalization are performed accurately. The choice of spatial resolution for track quantification studies is highly dependent on the amount of tracks the image is being sampled from. In our study, each subject was seeded with 5 million tracks, and we found that this led to reliable TDI maps with 1 mm3 isotropic voxels. Smaller voxel sizes led to maps with distributed empty voxels and greatly increased the computational cost during both TDI generation and normalization. Optimal settings for track quantification need to be standardized for it to become widely adopted.

This study was performed carefully and differences between groups were assessed conservatively. Our results demonstrate that TDI is affected in Parkinson's disease in specific areas of the brain, many of which have been previously linked to known dysfunctions. If these changes were global or randomly located they would not have been identified as significant by group-wise statistical tests.

Track-weighted imaging is inherently a non-local analysis technique. That is, biological changes or artifacts that are present in one region will impact the track-weighted measures sampled in other voxels traversed by the same tracks. This quality of track density imaging also likely contributed to our widely spatially distributed results. Even if the above issues are resolved, TWI may be more appropriate (and interpretable) in pathologies where there are only a small number of affected tracts.

The clinical relevance of neuroimaging findings when patient and control groups are compared is usually appreciated by testing for a simple relationship between clinical and imaging data. Our analysis did not identify any significant linear correlation between streamline density and individual performance in motor and cognitive tests: individuals with poorer performance did not show more severe imaging disturbances. This may suggest that the streamline density imaging method used in the present study is a biomarker of disease traits but may not be appropriate to monitor disease severity. Several other effects also need to be considered including compensatory mechanisms by exogeneous factors, such as anti-Parkinsonian medications at the time of clinical testing, endogeneous processes (i.e., neural network compensatory mechanisms), or a combination of both.

Potential future track quantification studies may want to consider using (i) outlier rejection methods like HOMOR (Pannek et al., 2012), (ii) lower maximum harmonic orders for CSD, or damped Richardson-Lucy deconvolution (Parker et al., 2013), (iii) templates created from subjects' FOD (Raffelt et al., 2011), (iv) anatomically constrained tractography (Smith et al., 2012) and filtering of tractograms given information from spherical-deconvolution (Smith et al., 2013), and (v) more reliable track-weighted imaging approaches (Besseling et al., 2012, Willats et al., 2014). Another method to consider may be direct comparison of FOD functions using apparent fiber density, which has proven successful in motor neuron disease (Raffelt et al., 2012). Global tractography approaches are also becoming feasible for widespread use (Reisert et al., 2010) and can provide more biologically plausible tractograms.

Conclusion

We have performed track density imaging within a large population of early-stage patients with Parkinson's disease and a group of matched healthy controls. Our cross-sectional analysis revealed increased streamline counts in Parkinson's disease within various white matter pathways traditionally involved in brainstem and extrapyramidal functions, limbic, and cognitive processes. The most remarkable result is the involvement of the nigrostriatal pathways extending dorsally from the ventrolateral part of the SN and terminating at the medial posterior putamen bilaterally.

Various sources of error and/or processing steps taken may have contributed to the increase in track density we observed in our patient population. At this point this confounds a biological interpretation of our results. What is clear is that track quantification can provide complementary information to standard diffusion analysis methods, and that it can be used to map both nigral and extranigral abnormalities in Parkinson's disease.