Sample size estimation
We calculated the optimal sample size of the corresponding expected effect size38 with G*Power 3.1.9.6 software. Based on our previous study15, we expected large effect sizes in the normalized volume of the SNc in the NM-MRI sequence comparing patients with controls (d = 1.20). Thus, the optimal sample size obtained for the expected effect size was 32 participants (with standard values α = 0.05, power = 0.8, one-tail tests).
Subjects
PD Patients fulfilled the UK Brain Bank Parkinson’s disease criteria39. In addition, patient selection was carried out considering their genetic status in order to enrich the genetically-determined PD groups and controls to match sex and age with them. HC individuals with no neurological disorders were recruited among their spouses. Written informed consent was obtained from all subjects, and the study was approved by the University of Navarra Research Ethics Committee.
Neurological and neuropsychological assessment
Patients underwent an interview covering demographic data, family history of neurological diseases, the Mini-Mental State Examination (MMSE)40, the Parkinson Disease-Cognitive Rating Scale battery (PD-CRS)41, the Unified Parkinson’s Disease Rating Scale (UPDRS-III)42, Modified Hoehn and Yahr Scale (HY)43, the Non-Motor Symptoms Scale (NMSS)44, and Edinburgh Handedness Inventory16. Levodopa equivalent daily dose (LEDD)45 was also calculated. Disease onset was established as the age when the Parkinsonian motor signs started self-reported when possible or reported by the caregiver. In addition, other clinical variables were assessed: “psychiatric symptoms” when depression or other psychiatric conditions were present, “hyposmia” (from NMSS item 28), and “other sleep disorders” (from NMSS items 3 and 5). RBD was recorded with the REM Sleep Behavior Disorder Single-Question Screen46. The UPDRS-III was assessed in an ON state right before the MRI scan.
Genetic analysis
Genomic DNA was isolated from leukocytes. Patients were screened for LRRK2 p.G2019S mutation by a custom-designed allele-specific PCR Taqman® assay and GBA gene analysis long-range PCR with confirmatory Sanger sequencing was performed47.
Neuroimaging analysis
MRI protocol and acquisition
Levodopa dosing was rescheduled in 22 patients to avoid dyskinesia during MRI scan. MRI scans were carried out on a 3-T MAGNETOM Skyra MRI scanner (Siemens Healthineers, Erlangen, Germany), using a 32-channel head array during a 32-min session for anatomical acquisition, 3D NM-sensitive sequence, and susceptibility-weighted images (SWI) datasets. The anatomical T1-weighted image was acquired with an MPRAGE sequence of 5 min. The following imaging parameters were employed: 1 mm-isotropic resolution, field of view (FOV) = 256 × 192 mm2, matrix = 256 × 192 voxels, 160 axial slices, repetition time (TR)/echo time (TE) = 1620/3.09 ms, Inversion Time = 650 ms, flip angle = 15º. The NM-sensitive sequences (i.e., NM-MRI) were obtained with a 3D-NM-sensitive T1-weighted turbo spin-echo sequence48 with the following parameters: repetition time/echo time, 34/4.91 ms; flip angle = 20º; 40 slices, 1 mm slice thickness, 0.2-mm gap, 512 × 408 acquisition matrix, 220 × 175 FOV, (voxel 0.6 × 0.6 × 1.0 mm), bandwidth 190 Hz/pixel, four averages. The SWI sequences (i.e., iron-MRI) were obtained by combining a long-TE high-resolution fully flow-compensated three-dimensional (3D) GRE sequence with filtered phase information in each voxel49 with the following parameters: TR/TE, 24/34 ms, flip angle = 10°, 22 cm FOV, 256 × 512 acquisition matrix, and 2 mm slice thickness (voxel-size 0.7 × 0.7 × 2.0 mm).
Image preprocessing
The slices were oriented orthogonally to the fourth ventricle floor and covered from the posterior commissure to the pons. Four excitations in the NM sequence and three excitations in the iron sequence were acquired and realigned offline to correct head movement. This image preprocessing was performed using SPM12 (Wellcome Trust Center for Neuroimaging, London, UK; http://www.fil.ion.ucl.ac.uk) and custom scripts in Matlab R2021a (Mathworks, MA, USA). In order to facilitate and standardize manual delineations, MR images were manually reoriented to adapt them to the orientation of the axial and mid-sagittal plane of a canonical T1 template image in SPM8.
Automatic image segmentation of brainstem structures
All the brainstem structures of interest were automatically segmented from the NM-MRI and iron-MRI sequences using the 3D-ABSP (Fig. 1). For this purpose, we created two static atlases of brainstem structures, one for the NM-MRI sequence and another one for the iron-MRI sequence, consisting of images of 32 HC and their corresponding manual annotations. The manual annotations were produced by an experienced neurologist (MAP), who delineated the structures of interest (SNc, LC, and brainstem in the NM-MRI sequences; whole iron deposit in SN, RN, and brainstem in the iron-MRI sequences) from all slices of the images of the 32 HC (Supplementary Fig. 8). Manual segmentation was blinded to the different groups, including controls.
The automatic segmentation of NM and iron-rich brainstem structures involved a multiresolution, three-step image registration of the target image with each of the atlas images, followed by a label fusion strategy to generate the final segmentation masks15. As quality-control of the segmentation, all brainstem structures of the 32 HC were segmented following a leave-one-out strategy (i.e., using an atlas elaborated from the remaining 31 HC), and the average Dice Similarity Coefficient (DSC) between the segmentation masks and their corresponding manual annotations for each brainstem structure was calculated as a quality-control segmentation accuracy score. After the quality of the segmentation was confirmed, the 71 PD images were segmented using the entire atlas composed of the images of the 32 HC.
Intra-subject spatial alignment of NM and iron sequences
Each subject’s individual NM and iron MRI images were aligned using a multiresolution rigid registration algorithm implemented with Elastix50. To assess the quality of the alignment, the DSC between the brainstem manual annotations of the original NM image and the transformed iron image was calculated.
Quantification of NM and iron in brainstem structures
The labels produced by the automatic segmentation of brainstem structures were used to quantify the amount of NM in SNc and LC, and the amount of iron in SNc and RN. Note that the intra-subject accurate alignment of both sequences allows applying the SNc segmentation mask obtained from the NM image to the iron image and constrains the iron measurement to the SNc. We measured the contrast ratio (CR), and the volume of each brainstem structure which was normalized to the total volume of gray matter (nVol). The gray matter volume was estimated with the get_totals function in Matlab. The CR was defined as the relative increase of the average intensity of the structure compared with the average intensity of the brainstem (normalized brightness). The automatically segmented structures were thresholded to measure the CR and nVol from hyperintense voxels in the case of NM15, and from hypointense voxels in the case of iron (Supplementary Figs. 9 and 10). CR values of iron were inverted in sign to facilitate interpretation (i.e., higher CR is translated into darker regions in the case of iron, whereas it means the region is brighter in the case of NM).
NM-iron interactions in the SNc were analyzed in-depth following an inter-subject image alignment strategy to create NM and iron spatial probabilistic maps and to display their specific spatial distribution patterns in the SNc. The topographical distribution of NM and iron content in this structure was thoroughly analyzed by fitting an elliptical-section cylinder which divided the SNc into four anatomical quadrants: Q1 corresponding to the ventral, Q2 mainly to the lateral, Q3 to the medial-rostral, and Q4 to the medial-caudal SNc (Fig. 1, Supplementary Fig. 11). All calculations were performed in the subject space.
Probabilistic maps of NM and iron content in brainstem nuclei
A reference HC image of the NM sequence with the highest DSC after the SNc segmentation was selected as the reference. Then the remaining NM sequence images were registered to the reference image following the same three-step registration framework described above.
Similarly, all the iron sequence images that were previously registered to their corresponding NM sequence image, underwent the same transformation to align them to the reference NM image. All brainstem structures’ segmentation masks underwent the same transformation as well. Hence, all the NM and iron sequence images were transformed into the same coordinate space, which allowed the creation of NM and iron content probabilistic maps using the segmentation masks of brainstem structures. Each voxel in the probabilistic map contained the sum of the segmentation masks of all the images (i.e., ‘1’s where the segmented structure of interest is present and ‘0’s where it is not) normalized to the number of subjects in the group for which the map was created (in the SNc N = 32 for HC and N = 63 for PD; in the LC N = 32 for HC, N = 39 for iPD and N = 24 for LRRK2-PD).
Spatial distribution of NM and iron in the SNc
A refined spatial analysis of the NM and iron content in the SNc was carried out by measuring the contrast ratio (CR) and normalized volume (nVol) in each anatomical quadrant and axial slice of the structure. To this end, NM and iron sequence images were aligned to the reference NM image as depicted in the previous section. Then, an elliptical-section cylinder was optimally fitted to the SNc, dividing the structure into four distinct anatomical quadrants. Q1 corresponded to the ventral SNc, Q2 mainly to the lateral SNc, Q3 to the medial-rostral SNc, and Q4 to the medial-caudal SNc. The automatic 3D segmentation was used to mask each quadrant in the elliptic cylinder, and CR and nVol per quadrant and slice were calculated for NM and iron signals (Supplementary Fig. 11).
Statistical analyses
Considering the distribution of the population studied, non-parametric (robust) statistical analyses were employed for group comparisons. One-way robust ANOVAs (WRS2 R package51; trimming value = 0.2) were run to examine group differences in the CR and nVol of brainstem structures in each sequence. Two-way robust mixed ANOVAs were used to evaluate the combined effects of group and region, group and side, and group and sequence. Post hoc t tests from the bootstrap version (number of samples = 500) were used. The effect size of post hoc comparisons was evaluated with a robust, heteroscedastic generalization of Cohen’s d52. The confidence intervals and the p values were adjusted for multiple testing53. Potential associations between brainstem MRI measures and demographic and clinical variables for the entire cohort of PD patients were evaluated using robust multiple regression analyses with MASS R package54. Robust multiple regression analyses with Huber M estimator (MASS R package) were applied to predict each of the brainstem MRI measures from demographic and clinical variables in PD subjects (i.e., age, sex, handedness, most affected side, years of education, disease duration, H&Y, UPDRS-III, LEDD, MMSE, PD-CRS, and presence of psychiatric disorders, hyposmia, RBD, or other sleep disturbances). To consider the amount and the interdependent nature of predictors, meaningful variables for each MRI parameter were obtained by means of feature selection with the Boruta algorithm55. Finally, robust regression analyses were repeated to predict each MRI brainstem measure from the meaningful variables suggested by the feature selection algorithm. Multiple-constraint hypotheses were conducted using robust Wald tests. Finally, the diagnostic performance of the 3D-ABSP was assessed by means of binary logistic regression analyses in Matlab. Bootstrapping with 1000 iterations was used to obtain AUC confidence intervals with 95% level of confidence. All tests were two-tailed, and p values <0.05 (corrected for multiple comparisons) were considered statistically significant.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.