EHS
EHS

# Cerebro-cerebellar motor networks in clinical subtypes of Parkinson’s disease

### Participants

One hundred-fifty-four patients with idiopathic PD were prospectively recruited at the Clinic of Neurology, Faculty of Medicine, University of Belgrade, Serbia, within the framework of an ongoing longitudinal project. Patients were included regardless of age, sex, education, family history of PD, age at disease onset, disease duration, side of onset, and levodopa equivalent daily dose (LEDD)42. Patients were excluded if they had: Hoehn and Yahr43 stage score >4; moderate/severe head tremor at rest; known PD-related genetic mutations (i.e., parkin and leucine-rich repeat kinase 2); dementia44; cerebrovascular disorders or intracranial masses on routine MRI; history of traumatic brain injury; any other major neurological and medical condition. Patients were assessed by clinical, cognitive/behavioral, and brain MRI evaluations (T1-weighted and resting-state functional MRI [rs-fMRI]) at baseline and clinically every year for a maximum of 4 years45,46. Eight PD patients were excluded from analysis due to incomplete MRI or movement artifacts.

An experienced neurologist blinded to MRI results performed clinical assessments45,46. Patients were examined in ON state (i.e., period when dopaminergic medication is working and symptoms are well controlled). Demographic, general clinical, and family data (sex, education, age, handedness, age at onset, side of onset, PD duration, and family history) were obtained using a semi-structured interview. LEDD was calculated. Disease severity was defined using the HY stage score and UPDRS. UPDRS was used to evaluate non-motor symptoms (UPDRS I), motor symptoms (UPDRS II), motor signs (UPDRS III) and motor complications (UPDRS IV). UPDRS III rigidity, axial, and bradykinesia subscores were also calculated. Severity of (FOG) was assessed using the FOG questionnaire. The presence of hallucinations was reported according to the UPDRS I subscore and of dyskinesia and fluctuations according to the UPDRS IV subscores. Cognitive and behavioral assessment was performed as previously described45,46. Non-motor symptoms other than cognitive impairment (gastrointestinal, urinary, olfactory, orthostatic and sexual dysfunctions) were assessed according to the Non-Motor Symptoms questionnaire (NMS-Q).

Sixty age- and sex-matched healthy controls without any neurological, psychiatric, or other disorders were also recruited among friends and relatives of patients and by word of mouth and performed the same study protocol.

Ethical standard committees on human experimentation of IRCCS Ospedale San Raffaele, Milan, Italy and Clinic of Neurology, Faculty of Medicine, University of Belgrade, Serbia approved the study protocol; all participants provided written informed consent prior to study inclusion.

### PD motor subtypes definition

PD patients were categorized in two motor subtypes (PD-TD and PD-PIGD) using the Movement Disorder Society Unified Parkinson’s Disease Rating Scale (MDS-UPDRS)47 as previously described2. Starting from the original sample (146 PD patients), 49 cases were classified as PD-TD and 83 as PD-PIGD, while 14 indeterminate patients were discarded from the analysis. We made the two PD groups more comparable to capture biologically relevant neuroimaging differences and avoid possible biases such as a disproportionate disease duration in PD-PIGD subjects. An “exact” matching analysis was applied in which age, sex, education, disease severity (UPDRS-III) and disease duration were used as covariates in PD group matching. Two groups of 32 PD-TD and 26 PD-PIGD patients were identified (Table 1 and Supplementary Tables 1 and 3). Analysis was performed using MatchIt package in R Statistical Software (version 4.0.3; R Foundation for Statistical Computing, Vienna, Austria).

Since PD-TD may convert into PD-PIGD along the disease progression, longitudinal clinical data of PD-TD cases were examined. PD-TD subjects were divided in those who converted to PIGD over 4-year follow-up (cPD-TD) and those who did not (ncPD-TD). PD-TD subjects presenting inconsistences or missing data in the clinical data reports during follow-up were excluded from this subanalysis.

### MRI analysis

#### Mapping motor activation in the cerebellum

Twenty-three right-handed independent healthy controls were recruited at the San Raffaele Scientific Institute by word of mouth and performed neuropsychological and MRI assessments (Supplementary Table 4). Subjects were excluded if they had: medical illnesses or substance abuse that could interfere with cognition; any major systemic, psychiatric, neurological, visual, and musculoskeletal disturbances; contraindications to undergo MRI examination; brain damage at routine MRI, including lacunae and extensive cerebrovascular disorders.

Healthy participants underwent a neuropsychological assessment. A blinded and experienced neuropsychologist performed a comprehensive cognitive evaluation including: Mini-Mental State Examination (MMSE)48; Digit span forward49 and Rey Auditory Verbal Learning Test (RAVLT) immediate and delayed recall50 and recognition51; Modified Card Sorting Test52; phonemic and semantic verbal fluency tests53; Attentive Matrices Test54; Trail Making Test (TMT)55 and digit span backward56; Freehand copying of drawings with and without landmarks50; Token Test57; Beck depression inventory (BDI)58; Apathy Rating Scale59; Snaith-Hamilton Pleasure Scale (SHAPS)60.

Brain MR images were acquired using a 3.0 Tesla Philips Intera scanner (Ingenia CX, Philips Medical Systems, Best, The Netherlands). The following sequences were acquired in all subjects: (i) three-dimensional (3D) T1-weighted sequence (repetition time (TR) = 7.1 ms, echo time (TE) = 3.2 ms, flip angle = 9°, 204 contiguous sagittal sections, thickness = 1 mm, field of view (FOV) = 256 × 240 mm, matrix = 256 × 240, voxel reconstruction 1 1 × 1 mm); (ii) 3DT2-weighted sequence (TR = 2500 ms, TE = 330 ms, flip angle = 90°, 192 contiguous sagittal sections, thickness = 1 mm, FOV = 256 × 256 mm, matrix = 256 × 258, voxel reconstruction 0.9 × 0.9 × 1 mm); and (iii) fMRI using a T2* weighted echo planar imaging (EPI) sequence (TE = 35 ms, TR = 1572 ms, flip angle = 70°, FOV = 240 × 240 mm, matrix = 96 × 94, 48 contiguous axial sections, thickness = 3 mm, acquisition time = 3 min and 57 s, voxel reconstruction 2.5 × 2.5 × 3 mm). fMRI scans were acquired during the performance of a “motor-task”, which consisted in alternated self-paced dorsal and plantar flexion movements of the feet with eyes closed, with the knees supported by a soft wedge and flexed of about 30°. To standardize the amplitude of movements, subjects were asked to reach a fixed wood bar with their insteps so that the dorsal flexion was about 90°. A block design (ABAB) was used, in which the activation A (lasting about 19 s) corresponded to the dorsal and plantar flexion, while during the resting period B (lasting about 19 s) subjects were asked to maintain their eyes closed, and each period was repeated 6 times (total duration 3′57″). Subjects performed the task according to auditory stimuli (“go” and “stop”) at the beginning and at the end of the movement. A visual signal (green/red light), visible to the operators outside the MRI scanner, was projected in order to monitor the correct time of task execution. Subjects were trained to perform the task outside the scanner and were carefully monitored visually by an observer inside the scanner room during acquisition to ensure a correct performance. Tasks were performed well by all the subjects.

Cerebellums were isolated from 3DT1-weighted images using the SUIT toolbox in SPM12 software (http://www.fil.ion.ucl.ac.uk/spm, Wellcome Trust Center for Neuroimaging, London). Individual white and gray matter segmented maps of cerebellum were then normalized into the SUIT atlas template using Dartel space. The normalization process uses the tissue segmentation maps, the white and gray matter segmentation maps resulting from the previous step of isolation of the cerebellum. The flowfield and affine transformations resulting from the normalization process were then used to bring different images from the anatomical space into SUIT space.

Task-based fMRI data were analyzed using SPM12 software. For each participant, all the volumes of fMRI data were realigned to the first one to correct for head movement (all study participants showed maximal head movements lower than 3 mm in each direction). Realigned functional data were then coregistered to the anatomical T1 images. All first-level analyses were performed at single-subject level. Signal variations of the BOLD effect associated with the execution of motor-task were evaluated voxel by voxel running the first-level General Linear Model (GLM). Specific effects were tested applying appropriate linear contrasts. Significant hemodynamic changes for each contrast were evaluated using Statistical Parametric Maps “t” (SPMt). Contrast images, resulting from first-level analysis, were then resliced into SUIT space applying flowfield and affine transformations obtained from normalization step. fMRI data in the SUIT space were then smoothed applying a 4-mm 3D-Gaussian filter. Second-level analysis was performed on smoothed fMRI data. One-sample t-test in SPM12 was performed to evaluate significant mean brain activations of healthy controls. A segmented cerebellum mask of healthy subjects was used as inclusion mask. All findings are shown at p < 0.001 uncorrected at the voxel level but only clusters passing a small volume correction for multiple comparisons (10 mm radius), cut-off value for significance p < 0.05, were presented. Cerebellar activations, resulting from second-level analysis, resulted into a 2 × 2 × 2 mm3 space. Finally, activations were normalized to 5 × 5 × 5 mm3 space, in order to obtain the motor cerebellar seed regions for the SFC analysis (Fig. 1A).

### MRI acquisition

Brain MRI scans were acquired on the same 1.5 Tesla Philips Medical System Achieva machine. Subjects were scanned between 10 and 11 a.m., i.e., treated PD patients were 90–120 min after their regular morning dopaminergic therapy administration (ON state). The following MR sequences were obtained: (i) dual-echo (DE) turbo spin-echo (SE) (TR = 3125 ms, TEs = 20/100 ms, echo train length [ETL] = 6, 44 axial slices, thickness = 3.0 mm, matrix size = 256 × 247, FOV = 240 × 232 mm2; voxel size, 0.94 × 0.94 × 3 mm, in-plane sensitivity encoding [SENSE] parallel reduction factor, 1.5); (ii) 3D sagittal T1-weighted Turbo Field Echo (TFE) (frequency direction = anterior-posterior, TR = 7.1 ms, TE = 3.3 ms, inversion time = 1000 ms, flip angle = 8°, matrix size = 256 × 256 × 180 [inferior-superior, anterior-posterior], FOV = 256 × 256 mm2, section thickness = 1 mm; voxel size = 1 × 1 × 1 mm, out-of-plane SENSE parallel reduction factor = 1.5, sagittal orientation); and (iii) gradient-echo (GRE) EPI for rs-fMRI (TR = 3000 ms, TE = 35 ms, flip angle = 90°, matrix size = 128 × 128, FOV = 240 × 240 mm2; slice thickness = 4 mm, 200 sets of 30 contiguous axial slices). During RS fMRI scanning, subjects were instructed to remain motionless, to keep their eyes closed, and not to think about anything in particular.

### Image preprocessing

rs-fMRI data processing was performed with the Data Processing Assistant for Resting-State toolbox (DPARSFA, http://rfmri.org/DPARSF61), based on Statistical Parametric Mapping (SPM12, http://www.fil.ion.ucl.ac.uk/spm), and the rs-fMRI Data Analysis Toolkit [http://www.restfmri.net62]. Preprocessing included the following: (1) removal of first four volumes of each raw rs-fMRI dataset to allow for T1 equilibration; (2) slice timing correction for interleaved acquisitions (the middle slice was used as the reference point); (3) head motion correction using a six-parameter (rigid body) linear transformation with a two-pass procedure (registered to the first image and then registered to the mean of the images after the first realignment); (4) spatial normalization to the Montreal Neurological Institute (MNI) atlas template with voxel size was set at 5 × 5 × 5 mm3 for computational efficiency; (5) removal of spurious variance through linear regression: including 24 parameters from the head motion correction step [6 head motion parameters, 6 head motion parameters one time point before, and the 12 corresponding squared items63], scrubbing with regression [signal spike regression as well as 1 back and 2 forward neighbors64] at time points with a frame-wise displacement (FD) > 0.5 mm65, linear and quadratic trends, global signal, white matter signal, and the cerebrospinal fluid signal; (6) spatial smoothing with a 4 mm FWHM Gaussian Kernel; and (7) band-pass temporal filtering (0.01–0.08 Hz) to reduce the effect of low frequency drift and high frequency noise66,67.

No participant had more than 2 mm/degree of movement in any of the six directions, and no more than 90 volumes removed during scrubbing (1/3 of the total volumes), ensuring at least 5 min and 30 s of functional data per individual.

### Stepwise functional connectivity analysis

SFC analysis is a graph theory metric that detects both direct and indirect functional couplings between a predefined seed region and other brain regions (Fig. 1B)16,68,69. SFC analysis aims to characterize regions connected to specific seed brain areas at different levels of link-step distances16,68,69. With such a framework, a step refers to the number of links (edges) that belong to a path connecting a node to the seed (or target) area. In SFC analysis, the degree of stepwise connectivity of a voxel j for a given step distance l and a seed area i ($$A_{ji}^l$$) is computed from the count of all paths that (1) connect voxel j and any voxel in seed area i, and (2) have an exact length of l. Each SFC matrix Al of size m-by-m can be recursively represented as follows:

$$A_l\left( {i,j} \right) = \left\{ {\begin{array}{*{20}{c}} {A\left( {i,j} \right)\left[ {i \ne j,l = 1} \right]} \\ {\mathop {\sum}\limits_{k = 1}^m {\left( {\frac{{A_{l – 1}\left( {i,k} \right) – \min \left( {A_{l – 1}} \right)}}{{\max \left( {A_{l – 1}} \right) – \min \left( {A_{l – 1}} \right)}}} \right)\left( {\frac{{A\left( {k,j} \right) – \min \left( A \right)}}{{\max \left( A \right) – \min \left( A \right)}}} \right)\left[ {i \,\ne\, j,l \ge 2} \right]} } \end{array}} \right.$$

(1)

Here, Al is the functional connectivity matrix with a step distance of l, and A is the correlation matrix after Fisher transformation. Matrices were then normalized between 0 and 1, keeping the final distribution of values intact while making them comparable across step distances. In this sense, a larger SFC degree under the step distance l indicates stronger paths connecting two voxels via link one, while a smaller degree indicates weaker connectivity paths. It is easy to see that SFC is also related to transition probabilities in an information flow analysis between the seed area and voxels across the rest of the brain. Given the lack of directionality information provided by rs-fMRI data, in SFC we did not include any restrictions about recurrent pathways crossing the seed regions multiple times. Therefore, we explored a wide range of link-step distances, from 1 to 20, to characterize the progression of the derived maps. However, as we found that the SFC patterns are topographically dissimilar between consecutive maps from steps one to three, and become stable for link-step distances above four, we only included maps from one to four steps in our results. SFC was applied to characterize pattern of rs-fMRI connectivity between the motor cerebellar seed region and the rest of the brain in PD-TD, PD-PIGD and matched 60 healthy controls (Fig. 1B).

### Network construction

Association matrices for each participant were computed by calculating the Pearson correlation between each voxel time course and every other voxel time course within a mask covering cortical and subcortical gray matter. To perform this analysis, preprocessed resting state images of each participant were previously converted to an N-by-M matrix, where N was the image voxels in MNI space, and M was the 200 acquisition time points. From this step, a 11,705 × 11,705 matrix of Pearson correlation coefficients (r-values) was obtained for each individual. Fisher z transformation was applied to r-values. Then, all negative correlations and positive correlations that did not reach a false discovery rate (FDR) correction67 of p-value < 0.05 were excluded from further analyses. Final association matrix included only significant positive associations, as positive connectivity has been proved to drive functional connectivity network topology in the human brain.

### Statistical analyses

Demographic, clinical, and cognitive/behavioral data were compared between groups using ANOVA models (for continuous variables) or Chi-square test (for all categorical variables). Two-sided p-value < 0.05 was considered for statistical significance. P-values were adjusted for multiple comparisons (R Statistical Software).

In order to identify regions demonstrating between-group SFC differences, voxel-wise analyses were performed using general linear models as implemented in SPM12 (Wellcome Department of Imaging Neuroscience, London, England; www.fil.ion.ucl.ac.uk/spm). Whole-brain two-sample t-test comparisons between groups were performed, including age and sex as covariates. A threshold-free cluster enhancement method combined with nonparametric permutation testing (5000 permutations) as implemented in the Computational Anatomy Toolbox 12 (CAT12, http://www.neuro.uni-jena.de/cat/) was used to detect statistically significant differences at p-value < 0.05, family-wise error (FWE) corrected.

The discriminatory power of clinical data and/or SFC at different link-step distances in distinguishing cPD-TD from ncPD-TD patients was assessed using ROC curve analyses (Fig. 1C). Statistical significance level was set at p-value < 0.05. The area under the curve (AUC), as derived measure of accuracy, was calculated (R Statistical Software).

EHS