Study samples
Network metrics were computed for previously validated motor and cognition-related PD topographies described elsewhere1,3. Group differences in the measures, as well as the effects of disease progression, genotype, and treatment were explored in previously published datasets referenced below (Table S1A-S2E). PD patients in each sample were diagnosed according to UK Parkinson Disease Society Brain Bank criteria35. The subjects were evaluated at the time of imaging according to the Unified Parkinson’s Disease Rating Scale (UPDRS)36 at least 12 hours after the last medication dose. Ethical permission for these studies was obtained from the Institutional Review Boards of the participating institutions. Written consent was obtained from each subject following detailed explanation of the procedures.
Study cohorts
1. Comparison of network metrics in PD and healthy control subjects
(a) [18F]-fluorodeoxyglucose (FDG PET). We studied 96 PD patients and 22 age- and gender-matched healthy control (HC) subjects (HC1, n=22) who were scanned with FDG PET at North Shore University Hospital as part of long-term diagnostic project15 (Table S1A). For validation, we studied an independent cohort of 146 PD and 39 HC subjects who were scanned with FDG PET as part of two subsequently published ascertainment studies16,17 (Table S1B).
(b) rs-fMRI. To further substantiate the finding in single case measurements of PDRP assortativity, we studied a cohort of 20 PD and 20 age- and gender-matched HC subjects (HC2, n=20) at North Shore, who were scanned with resting-state fMRI (rs-fMRI)7 (Table S1C). To validate the findings, we additionally analyzed baseline rs-fMRI scans from 102 PD and 18 HC who were studied as part of Parkinson’s Progression Markers Initiative (PPMI; https://www.ppmi-info.org/)37 (Table S1C).
2. Change in network metrics with disease progression
(a) Longitudinal progression. We studied a longitudinal cohort of early PD patients (n=15) who were scanned with FDG PET at baseline, 24 months and 48 months, and a matched group of healthy control subjects (HC3; n=15) (see Table S2A). Limited metabolic data from these subjects have appeared previously18.
(b) Comparison of groups with increasing symptom duration. We compared three age- and gender-matched groups of PD patients with motor symptoms of short (0-4 years; n=20), intermediate (5-9 years; n=20), and long (≥10 years; n=20) duration, and age- and gender-matched healthy control subjects (HC4; n=20) (see Table S2B).
3. Effects of PD genotype on network metrics
We studied patients with PD associated with either the LRRK2-G2019S mutation (PD-LRRK2; n=14) or GBA1 variants (PD-GBA; n=12), and matched healthy control subjects (HC5; n=14) (Table S2C). Clinical details and limited network data from these patients have appeared previously8.
4. Effects of levodopa infusion on network metrics
We studied PD subjects (n=14) who underwent metabolic imaging in the unmedicated state (off) and during intravenous levodopa infusion (on), and a matched group of healthy control subjects (HC6; n=14) (Table S2D). In each patient, the levodopa dose was titrated to achieve maximal clinical benefit, without eliciting involuntary movements. Experimental details and limited metabolic data from these patients have appeared previously23.
5. Effects of subthalamic gene therapy and sham surgery on network metrics
We analyzed baseline and 12-month FDG PET data from PD patients who were randomized to either subthalamic AAV2-GAD (n=16) or sham surgery (n=21) as part of a double blind surgical trial, and a matched group of healthy control subjects (HC7; n=22)24 (Table S2E). The results of graph analysis of the gene therapy-induced network topography have appeared previously4.
Imaging
FDG PET imaging
Study participants at North Shore underwent metabolic imaging with FDG PET on a GE Advance tomograph (General Electric, Milwaukee, WI) at The Feinstein Institutes for Medical Research (Manhasset, NY) as described in detail elsewhere6,38. Study participants from other centers underwent metabolic imaging on their respective FDG PET platforms16,17,24. In all centers, subjects fasted overnight prior to PET imaging; anti-parkinsonian medications were withheld for at least 12 hours before the scanning session. Scans from each subject were realigned and spatially normalized to a standard Montreal Neurological Institute (MNI)-based PET template and smoothed with an isotropic Gaussian kernel (10 mm) in all directions to improve the signal-to-noise ratio. Image processing was performed using Statistical Parametric Mapping (SPM5) software (Wellcome Trust Centre for Neuroimaging, Institute of Neurology).
MR imaging
All participants received MR scans including structural and functional imaging. During the resting-state sequence they were instructed to rest quietly with eyes open and not to sleep.
Study participants at North Shore were scanned on a General Electric 3.0 Tesla Signa HDxt scanner. For each subject, a 3D T1 image was acquired using a SPGR protocol with repetition time between 2 RF pulses (TR) = 7.6 ms, echo time (TE) = 2.98 ms, TI = 650ms, field of view (FOV) = 240 mm, flip angle (FA) = 8°, and voxel size = 1×1×1 mm3. For rs-fMRI, 240 volumes × 40 slices were acquired for 8 minutes with TR = 2000 ms, TE = 28 ms, FOV = 240 mm, FA = 77° and voxel size = 3 mm. Further details on scanning parameter can be found elsewhere7.
PPMI study participants were scanned on a 3T Siemens Trio Tim scanner. For each subject, a 3D T1 image was acquired using a MPRAGE protocol with repetition time (TR) = 2300 ms, echo time (TE) = 2.98 ms, field of view (FOV) = 256 mm, flip angle (FA) = 9°, and voxel size = 1×1×1 mm3. For rs-fMRI, 210 volumes × 40 slices were acquired for 8 minutes and 24 seconds with TR = 2400 ms, TE = 25 ms, FOV = 224 mm, FA = 80° and voxel size = 3.3 mm. Further details on scanning parameter can be found on the website (http://www.ppmi-info.org/wp-content/uploads/2017/06/PPMI-MRI-Operations-Manual-V7.pdf). Image preprocessing was performed using the FMRIB Software Library (FSL; www.fmrib.ox.ac.uk/fsl) for rs-fMRI scans from both cohorts. It included motion correction, brain extraction, spatial smoothing (kernel = 8 mm; FWHM) and temporal high-pass filtering (cutoff frequency = 1/100 Hz). Rs-fMRI volumes were registered to the individual patient’s structural T1 image and then to the standard MNI 152 template. Further details of image preprocessing can be found elsewhere7.
Network analysis
1. PD-related networks
We used previously validated PD-related metabolic patterns associated with motor and cognitive dysfunction, known respectively as PDRP and PDCP3,33.
These patterns were identified using spatial covariance analysis of FDG PET data from PD patients and HC subjects as described previously6,31,39. Based on prior graphical analysis, we used nodal valence (region weight sign) to partition both networks into core and periphery subgraphs22. For the PDRP (38 nodes according to the AAL atlas40 (see below; Table S3A), the 20 nodes with positive metabolic valence (i.e., those with region weights ≥ 1.0) were found to form a discrete topological core zone8,22, whereas the remaining 18 nodes had negative valence (i.e., region weights ≤ -1.0) and comprised the network periphery8. For PDCP (35 nodes; Table S3C), the core was composed of the 16 nodes with negative valence8, which corresponded to the cortical core of the default mode network (DMN)3,31,41. The assignment of each PDRP and PDCP node to the core or periphery of the respective network is given in Tables S3A and S3C. For both networks, assignments based on metabolic valence agreed well with the results of modularity maximization and information theoretic community detection algorithms8.
For the rs-fMRI datasets, we used the fPDRP, a PD-related motor pattern topographically similar to the PET-based PDRP7. The fPDRP was identified using independent component analysis of rs-fMRI data from patients and control subjects in conjunction with bootstrap resampling7. It is composed of 39 nodes according to the AAL atlas40 (Table S3B).
2. AAV2-GAD-related pattern (GADRP)
We utilized the previously reported AAV2-GAD-related metabolic pattern (GADRP) identified in trial participants scanned with FDG PET at baseline and after subthalamic gene therapy4,24. This network was extracted from the metabolic data using ordinal trends canonical variates analysis (OrT/CVA), a form of supervised principal component analysis42. The details of this analysis and its results are provided elsewhere4. The GADRP was composed of 14 nodes according to the AAL atlas (Table S3D). Because of its relatively small size, this network could not be reliably partitioned; graphical analysis was therefore performed on the network as a whole. Preliminary graph theory analysis of the GADRP has appeared previously4. In the current study, we extended the analysis to the PDRP and additionally measured graphical metrics for this space in the various groups and time points.
Graph theory and network metrics
Defining nodes and edges
As noted above, for each network, nodes were defined based on the AAL atlas40, in which we parcellated the 3D image of the whole brain, normalized to MNI space. This produced 95 standardized anatomical regions-of-interest (ROIs) as described previously4,8,22. For the significant clusters identified by voxel-wise network analysis, we identified corresponding AAL nodes (see Table S3A-D for listings). For each node, we computed normalized metabolic activity for FDG PET scans and mean time series for rs-fMRI scans.
Different network spaces were used in the various analyses. Most of the studies were analyzed in PDRP space for FDG PET6 and fPDRP space for rs-fMRI data7. As mentioned above, comparison of PD genotypes (PD-LRRK2 vs. PD-GBA) relied on prespecified subnetworks of the PDRP and PDCP (i.e., network core and periphery) defined in an earlier study8. Lastly, in addition to the PDRP, we used the treatment-induced AAV2-GAD gene therapy network (GADRP) reported previously4 to assess the changes in graph metrics that occurred with this intervention and with sham surgery in this network space and in the PDRP. The details of each of these networks are provided in Supplementary Tables S3A-S3D.
In this study, we used AAL ROI data from FDG PET (globally normalized regional metabolic activity) and rs-fMRI (filtered regional mean time series) to construct matrices of pairwise correlations at the group and subject levels. For group-level analysis of the FDG PET scans, we used bootstrap methods (in-house Matlab script; MATLAB R2020a) to generate 100 samples for each group. For each iteration, we computed pairwise nodal correlation coefficients (Pearson correlations). The median values of the iterates (100 bootstrap correlation estimates) were used to create an adjacency matrix for the network/subgraph in each group. For subject-level analysis of the rs-fMRI scans, regional mean time series data were filtered using a continuous Gabor wavelet transform in the 0.01-0.12 Hz frequency range (cwt function in MATLAB R2020a). We then computed pairwise nodal correlation coefficients to create an adjacency matrix for the network in each subject/group. These calculations were performed using the Wavelet Toolbox and Statistics and Machine Learning Toolbox in MATLAB R2020a.
Network metrics
To assess group differences in connectivity parameters within the relevant network spaces, we used undirected graphs for simplicity in hypothesis testing and computational ease. The following metrics were computed on weighted graphical links:
- Assortativity coefficient: the correlation coefficient between the degrees of all nodes on two opposite ends of a link11–13. The details of calculating the assortativity coefficient are described elsewhere13. For a given network, we consider assortativity in a group of subjects to be increased if the assortativity coefficient, which is described as assortative for positive values, neutral for ≈ 0, or disassortative for negative values, is significantly elevated compared to values for the same network computed in a different group. Analogously, assortativity is considered to be reduced when the coefficient computed for the same network is lower than in a comparison group.
- Degree centrality: the number of connections (edges) within the network or subgraph, divided by the total number of nodes in the same space.
- Clustering coefficient: a measure of the likelihood that the nearest neighbors of a node will also be connected.
- Characteristic path length: the shortest path length between two nodes averaged over all pairs of nodes in a given network. High characteristic path length implies less efficient information transfer through the network43,44.
- Small-worldness: the ratio of clustering coefficient to characteristic path length, normalized to corresponding parameters from an equivalent random graph41. This measure quantifies the ratio of segregation to integration of information sources in the network space.
These parameters were computed using the Brain Connectivity Toolbox44 and an in-house Matlab script (MATLAB R2020a). We present assortativity and the other network metrics over a range of connectivity thresholds as described previously4,8. In the current study, thresholds ranged from r=0.3 to 0.6, at 0.05 increments, corresponding to graph densities between 25% and 60%4,8. At lower thresholds (r<0.3, graph density >60%), however, group differences may be difficult to discern because of the inclusion of random, non-specific links. At higher thresholds (r>0.6, graph density <25%), graphs may disconnect and distort group comparisons. By plotting the results over a range of thresholds, we demonstrated that group differences in a given metric were robust beyond two or three adjacent levels.
Stability of network metrics
To evaluate test-retest reliability for PDRP network metrics, the 8-min rs-fMRI acquisitions in the North Shore PD cohort (n=20) were divided into two 4-min scans, which were randomly assigned to the test or retest conditions as described elsewhere7. Test-retest reliability was evaluated by computing the intraclass correlation coefficient (ICC) and 95% confidence interval (CI) across subjects and correlation thresholds.
Visualization
To visualize assortativity effects along a graph, we generated 2D displays of the joint probability distribution of remaining degree (i.e., the total degree centrality for each node, minus the connected edge) for random pairs of nodes in the network space11,13. This was done using the bivariate Gaussian copula with the correlation matrix S=[1 ; 1]. Poisson marginal distributions ( =1000 random numbers, mean=λ) were computed for each variable, where is the degree-degree correlation coefficient and λ is the average degree of the network11,13. These two parameters were estimated from the empiric data from each group and network. The simulation steps were performed using an in-house Matlab script (Statistics and Machine Learning Toolbox, MATLAB R2020a).
To visualize differences in connectivity patterns across groups or conditions, graphs were displayed at the minimum threshold (Level 1, r=0.30) for a given network. For assortative configurations (mean ρ>0), an exemplar with assortativity coefficient between the mean value and the mean +0.5 sd was selected from the bootstrap samples described above. For disassortative configurations (mean ρ<0). an exemplar was selected with assortativity coefficient between the mean value and the mean –0.5 sd. Graphical configurations were constructed using an in-house Matlab script (mathematics/graph and network algorithms toolbox, MATLAB R2020a).
Lastly, to assess differences in nodal organization between groups and/or conditions, we ordered the nodes hierarchically according to degree centrality. For each configuration, we computed the Spearman rank-order correlation coefficient for the median graph at minimum threshold (r=0.30). Significant correlations (p≤0.05) reflected similar nodal hierarchy across groups or conditions. Uncorrelated rank order implied the presence of structural differences in the two graphs.
Statistical analysis
Network expression values were compared across groups using student’s t-test or one-way analysis of variance (ANOVA) with Tukey-Kramer HSD post-hoc tests for multiple comparison correction. For the graph analysis, the bootstrapped data or subject data were used to assess group differences in the network parameters for the relevant subgraphs. Group differences in each of the parameters were evaluated using the general linear model with post-hoc Bonferroni tests across graph thresholds. These analyses were performed using IBM SPSS Statistics for Windows, version 21 (IBM Corp., Armonk, N.Y., USA) or GraphPad Software (Version 7.0, La Jolla California, USA). Results were considered significant for p<0.05, corrected.
Data availability
Deidentified data will be made available on reasonable request from interested investigators for the purpose of replicating results.