A Radiomics Approach to Predicting Parkinson’s Disease by Incorporating Whole-Brain Functional Activity and Gray Matter Structure

Parkinson’s disease (PD) is a progressive, chronic, and neurodegenerative disorder that is primarily diagnosed by clinical examinations and magnetic resonance imaging (MRI). In this study, we proposed a machine learning based radiomics method to predict PD. Fifty healthy controls (HC) along with 70 PD patients underwent resting-state magnetic resonance imaging (rs-fMRI). For all subjects, we extracted five types of 6664 features, including mean amplitude of low-frequency fluctuation (mALFF), mean regional homogeneity (mReHo), resting-state functional connectivity (RSFC), voxel-mirrored homotopic connectivity (VMHC) and gray matter (GM) volume. After conducting dimension reduction utilizing Least absolute shrinkage and selection operator (LASSO), fifty-three radiomic features including 46 RSFCs, 1 mALFF, 3 mReHos, 1 VMHC, 2 GM volumes and 1 clinical factor were retained. The selected features also indicated the most discriminative regions for PD. We further conducted model fitting procedure for classifying subjects in the training set employing random forest and support volume machine (SVM) to evaluate the performance of the two methods. After cross-validation, both methods achieved 100% accuracy and area under curve (AUC) for distinguishing between PD and HC in the training set. In the testing set, SVM performed better than random forest with the accuracy, true positive rate (TPR) and AUC being 85%, 1 and 0.97, respectively. These findings demonstrate the radiomics technique has the potential to support radiological diagnosis and to achieve high classification accuracy for clinical diagnostic systems for patients with PD.


INTRODUCTION
Parkinson's disease (PD) is a major neurodegenerative disease influenced by both genetic and environmental factors (Halliday et al., 2014). As the second most common neurodegenerative disorder, PD is characterized by the degeneration of dopamine-producing cells in the brain resulting in motor symptoms and nonmotor features (Mhyre et al., 2012). Available diagnostic tools are better at detecting motor symptoms than nonmotor symptoms. The neural and pathophysiologic mechanisms to predict the progression of PD remain unclear and discovering the psychobiological markers is the key research priority. Understanding the inner working mechanisms of PD is one of the most intriguing scientific questions. Studies in neuroscience strongly suggest intervention during early therapeutic windows (Vu et al., 2012;Tibar et al., 2018). Although positron emitted topography/computed tomography is accurate (Meles et al., 2017) the diagnosis of PD at present is mainly dependent on clinical features and scores.
In recent years, neuroimaging has been increasingly employed to aid the early diagnosis of PD. A variety of neuroimaging technologies including functional magnetic resonance imaging (fMRI), structure MRI (sMRI), positron emission tomography (PET) and electroencephalography (EEG) have been widely adopted. Among these, resting-state functional MR imaging (rs-fMRI) is regarded as a promising technique for precisely locating the abnormal spontaneous activities in neuropsychological disease . Several rs-fMRI based methods including regional homogeneity (ReHo), amplitude of low frequency fluctuations (ALFF), and functional connectivity (FC) provide a task-free approach to explore spontaneous brain activity and connectivity among networks in different brain regions of PD patients. Application of these techniques provides new insights in prediction, early diagnosis and differential diagnosis of PD. Previous rs-fMRI studies using ALFF found specific frequency band of ALFF for PD  and detected significant alterations of ALFF in the prefrontal cortex and subcortical regions in PD patients (Xiang et al., 2016). Frequency domain analyses of ALFF revealed decreased ALFF in the putamen, parieto-temporo-occipital cortex, thalamus, cerebellum, and several occipital regions, while increased ALFF values were detected in the caudate and several temporal regions . A large sample study of 109 PD patients found distinguishing frequency bands and neural modulations in the brainstem and striatum correlated with the dose of levodopa and bradykinesia subscale scores (Hou et al., 2014).
Along with these distinctive changes of ALFF, spontaneous ReHo analysis in rs-fMRI studies also achieved considerable progress in examining early onset and late-onset PD (Wu et al., 2009;Yang et al., 2013;Sheng et al., 2016). ReHo alteration in the early phase of PD showed a low level of local coherence in the right primary sensory and positive correlation with disease duration (Choe et al., 2013). A 2-year longitudinal PD study of multimodal MRI using ReHo and voxel-based-morphometry (VBM) observed a progressive decrease of ReHo values in the sensorimotor cortex, default-mode network (DMN), and the left cerebellum, but increased ReHo in the supplementary motor area (SMA), bilateral temporal gyrus, and hippocampus (Zeng et al., 2017). A meta-analysis using ALFF and ReHo found consistent decreased activity in the putamen for PD patients that could serve as an independent validation of rs-fMRI (Wang et al., 2018). Another ALFF and ReHo based study demonstrated the disturbed DMN, SMA, basal ganglia (BG), and posterior cerebellar lobule in cognitively normal PD as compared with healthy controls (Harrington et al., 2017). These multilevel characteristics of rs-fMRI could effectively improve the discrimination accuracy of diagnosis.
Although previous rs-fMRI studies revealed widespread abnormal intrinsic networks in line with the pathophysiology of PD, these findings and biomarkers have not been extensively used for diagnosis, prediction or prognosis of PD in daily clinical practice. In recent years, a method called radiomics that extracts large amount of features from radiographic medical images into high-dimensional mineable data using data-characterization algorithms has received considerable attention, particularly in clinical oncology diagnosis (Valladares et al., 2020). Radiomics analysis employs multimodality medical images and machine learning techniques to extract many quantitative characteristics as objective, sensitive biomarkers of disease stage to potentially detect treatment effects . Applications of radiomics approach in the neurodegenerative and mental disorder disclosed the heterogeneity characteristics with a high accuracy that facilitate individualized diagnosis in patients with Alzheimer's disease, autism spectrum disorder and schizophrenia (Hofmann-Apitius et al., 2015;Salvatore et al., 2019;Wang et al., 2019). The radiomics technology that integrated the advantages of various models has been utilized to extract the characteristics for automated diagnosis of early PD and quantifying PD severity. These methods consist of voxel-based method (VBM), diffusion tensor imaging (DTI), functional connectome and connectivity measures among others (Shinde et al., 2019). A radiomics analysis of longitudinal Single-photon Emission Computed Tomography (PSECT) images demonstrated radiomic features significantly increased the prediction accuracy and were proved to be effective prognostic biomarkers of PD (Rahmim et al., 2017). More recently, a radiomics of deep neural nets on neuromelaninsensitive MRI demonstrated a test accuracy of 85.7% and revealed the substantia nigra pars compacta abnormalities in PD discriminating from atypical PD (Shinde et al., 2019). Another radiomics study of quantitative susceptibility mapping were shown to assist the diagnosis of idiopathic PD (Cheng et al., 2019). A classifier for early PD with an accuracy of 86.96% was identified from SVM training by extracting characteristics including ALFF, ReHo and RSFC from the gray matter (GM), white matter (WM) and cerebrospinal fluid (CSF) (Long et al., 2012).
Considering the above-mentioned radiomics approaches in existing PD studies, we aimed to utilize both rs-fMRI and sMRI to extract radiomic features including whole-brain functional activity (i.e., ALFF and ReHo), connectivity (i.e., RSFC and VMHC) and gray matter (GM). Our goal was to discover more effective biomarkers and to eventually develop an automated classification framework of early diagnosis for PD patients.

Participates and Clinical Evaluation
This study was approved by the Medical Research Ethical Committee of Nanjing Brain Hospital (Nanjing, China) in accordance with the Declaration of Helsinki, and written informed consent was obtained from all subjects. Seventy PD patients and fifty healthy controls (HC) were recruited. All the demographic characteristics and clinical symptom ratings were collected before MRI scanning and all patients were in the ON state during the MRI scan.
All subjects underwent a complete neurological and psychological status assess, and a review of medical history records. Mini-mental state examination (MMSE) was used to evaluate cognition. The severity of depression was quantified using the Hamilton Depression Scale (HAMD). The neurocognitive tests were administered to each participant individually by a professional appraiser in the neuropsychological research center. All HC participants were interviewed to rule out the presence including current or past psychiatric illness, history of psychiatric illness in first-degree relatives and/or current or past significant medical or neurological illness.
The demographic and clinical data of patients with PD and HC were compared using a Fisher's exact test (for sex), analysis of variance (ANOVA) (for age, education, MMSE and HAMD). The level of significance was set at p < 0.05.

Image Data Acquisition
Image data were acquired using a Siemens 3.0-Tesla signal scanner (Siemens, Verio, Germany) in the department of radiology within Nanjing Brain Hospital. Functional imaging data were collected transversely by using a gradient-recalled echo-planar imaging (GRE-EPI) pulse sequence with the following configurations: TR/TE = 200 ms/30 ms, flip angle = 90 • , matrix = 64 × 64, FOV = 220 × 220 mm, thickness/gap = 3.5/0.6 mm, in-plane resolution = 3.4 × 3.4 mm, slices = 31. For each subject, a total of 140 volumes were obtained, resulting in a total scan time of 280 s. High resolution anatomical images were acquired using a T1 fluid attenuated inversion recovery (FLAIR) sequence (TR/TE = 2530/3.34 ms, flip angle = 7 • , matrix = 256 × 192, FOV = 256 × 256 mm, slice thickness/gap = 1.33/0.5 mm, 128 slices covered the whole brain). The subjects were instructed to keep their eyes closed, relax their minds and remain as motionless as possible during the data acquisition. Rubber earplugs were used to reduce noise, and foam cushioning was used to fix the head to reduce motion artifacts. The MR images were retrieved from the archive by two experienced neuroradiologists (QH and XW).

Data Preprocessing
Image preprocessing procedure was carried out Data Processing Assistant for Resting-State fMRI 1 based on Statistical Parametric Mapping (SPM12 2 ) operated on the Matlab platform. The following steps were applied to the image data. For each subject, we first discarded the first five time points for signal equilibrium when subject was still adapting to the scanning noise. The remaining 135 images underwent slice-timing correction using the middle slice as the reference frame and head motion correction by regressing out 6 head motion signals (displacement on x, y, and z direction and 3 angular motion). Four subjects with more than 2.5 mm maximum displacement in any of the three dimensions or 2.5 • of any angular motion were removed. Next, T1-weighted structural images of each subject were coregistrated to the resulting functional images followed by images being segmented into gray matter (GM), white matter, and cerebrospinal fluid using a new segment and DARTEL segmentation algorithm. The functional images were spatially normalized to the Montreal Neurological Institute (MNI) space with 3 × 3 × 3 mm cubic voxels. After normalization, images were spatially smoothed with a 4 mm full width at half maximum (FWHM) Gaussian kernel and detrended using linear, quadratic or higher order polynomial algorithms. Note that when calculating Regional Homogeneity (ReHo), the smoothing procedure was omitted for maintaining the measuring accuracy. After smoothing and detrending, we further regressed out nuisance covariates including the Friston 24 motion parameters (Friston et al., 1996), white matter, global signals and cerebrospinal fluid signals and applied temporal filter (0.01 Hz < f < 0.08 Hz) to diminish high-frequency noise.

ReHo Analysis
The method of Regional Homogeneity (ReHo) (Zang et al., 2004) was proposed to analyze characteristics of regional brain activity and to reflect the temporal homogeneity of neural activity. It has been pointed out that some preprocessing methods especially spatial smoothing R-fMRI time series may significantly change the ReHo magnitudes (Zuo et al., 2013). To get rid of this potential issue, preprocessed rs-fMRI data without the spatial smoothing step was used for calculating ReHo. All individual ReHo maps were computed and then spatially smoothed with a 4 mm FWHM Gaussian kernel. In particular, we focused on the mReHo maps obtained by dividing the mean ReHo of the whole brain within each voxel in the ReHo map. We further segmented the mReHo maps and extract all the 112 ROI signals based on the Harvard-Oxford atlas (HOA) using the Resting-State fMRI Data Analysis Toolkit, REST 3 .

ALFF and VHMC Extraction
Slow fluctuations in activity are fundamental features of the resting brain for determining correlated activity between brain regions and resting state networks. The relative magnitude of these fluctuations can discriminate between brain regions and subjects. Amplitude of Low Frequency Fluctuations (ALFF) (Zang et al., 2007) are related measures that quantify the amplitude of these low frequency oscillations. Leveraging the preprocessed data within the frequency range between 0.01 and 0.1 Hz, we calculated individual ALFF maps and the mALFF maps by dividing the mean ALFF of the whole brain within each voxel in the ALFF maps. Using the HOA, we ended up with 112 mALFF values after extracting the ROI signals based on the mALFF maps.
Voxel-Mirrored Homotopic Connectivity (VMHC) quantifies functional homotopy by providing a voxel-wise measure of connectivity between hemispheres. VMHC calculates the connectivity between each voxel in one hemisphere and its mirrored counterpart in the other (Zuo et al., 2010). By segmenting the VMHC maps via HOA, we also got 112 VHMC values.

RSFC and GM Volume Extraction
Resting-state functional connectivity (rsFC) analysis is an effective method for estimating spontaneous functional activity and measuring the temporal correlation between spatially remote neurophysiological events. The preprocessed rs-fMRI images were segmented into 112 ROIs according to HOA. After averaging the rs-fMRI time courses of all the voxels within each ROI, the mean time series of each ROI were acquired. We performed Pearson's correlation analysis on each pair of ROI time series (i.e., 112 × 111/2 = 6216 pairs in total). The 6216 correlation coefficients were then transformed into z-scores by Fisher's z transformation and retained as the RSFC metrics. Based on the preprocessed structural images, we also extracted GM volumes of these 112 ROIs using the HOA as masks.

Feature Selection and Model Validation
Our candidate features consist of all the aforementioned metrics including ReHo, mALFF, VHMC, RSFC and GM volume along with all the clinical characteristics. To build our model, we first randomly split our dataset into training set and testing set while maintaining the PD:HC ratio, where 93 subjects were used as the primary cohort for feature selection and model training. The remaining 23 subjects were treated as validation cohort for examining the selected features. All steps of feature selection and model training were only based on and performed in the training dataset.
Note that our goal was to identify the most significant variables that could discriminate PD patients from healthy controls. However, as we had a total of 6669 features and a comparably much smaller sample size, the dimension reduction was necessary to improve the accuracy in the later step of building the machine learning model for classification . Hence, we first performed the nonparametric Mann-Whitney U test on each feature between the PD patients and the healthy controls and kept the variables with P-value larger than 0.1. Since the Mann-Whitney U test does not require the data to be normally distributed, we adopted the procedure as the first step to filter the features. Next, to avoid the possible presence of Simpson's Paradox caused by multicollinearity, where a predictor appears to be significant by itself, but this observation disappears or the direction reverses when other predictors are added. Therefore, whenever we spotted an absolute value of pairwise correlation between two features that was larger than 0.5, we removed the feature with larger average absolute correlation. Finally, to further reduce the burden of high dimensionality imposed on the model training, we used the least absolute shrinkage and selection operator (Lasso) procedure that assigns a penalty to the coefficients and eliminates variables with zero coefficient value. We used 10-fold cross validation to obtain the optimal penalty parameter for Lasso and retained the features with nonzero regression coefficients.
The methods of supper vector machine (SVM) and random forest were implemented for classifying subjects based on the selected features. Given the relatively small sample size, we used the linear kernel when fitting SVM, which was conducted by specifying the type parameter to be "linear" in the svm function. The tuning parameters for both methods were selected via 10-fold cross-validation by using the trainControl function with the method parameter being "repeatedcv" and the train function with the method parameter specified as either "cforest" or "svmLinear2" for random forest and SVM, respectively. The performance of these two different machine learning methods in the training and test sets were later compared and visualized according to different metrics including accuracy, true positive rate, false positive rate, receiver-operating characteristic (ROC) curve and area under curve (AUC). The complete statistical analysis was conducted in R 3.5.0. Specifically, packages "e1071, " "randomForest", "glmnet, " "caret" were employed for running the SVM, random forest, Lasso and for cross validation, respectively. The flowchart of this study is presented in Figure 1.

Clinical Characteristics
In Table 1, we provided the complete demographic and clinical information for all subjects participated in this study. No significant difference was observed with respect to the gender, age, education and MMSE score between PD patients and HCs, while significant difference was detected for HAMD between these two groups. In particular, for PD patients, the HAMD scores (11.0 ± 6.9) were significantly higher than those for HCs (2.1 ± 2.3).

Feature Selection
After the first step of Mann-Whitney U test, 6669 features containing metrics from rs-fMRI, sMRI and clinical information reduced to 993 features. Next, the procedure of excluding variables with absolute correlations larger than 0.5 removed 628 features with a total of 365 features remaining. Last, 54 features including (46 RSFCs, HAMD, 1 mALFF, 3 mReHos, 1 VMHC and 2 GM volumes) with nonzero coefficients obtained from the logistic regression with Lasso penalty were retained as the final metric set to be used for binary classification. In Table 2, we listed these 46 RSFCs and the respective connected brain regions indexed in the HOA template. The related brain regions of RSFCs were primarily located in the executive control network (ECN), default mode network (DMN), affective network (AN), visual network (VIN) and sensorimotor network (SMN) (Figure 2). The seven features were mALFF of the left superior temporal gyrus, posterior division, mReHo of the left parahippocampal gyrus, posterior division, mReHo of the right thalamus and left pallidum, VMHC of the right temporal fusiform cortex, anterior division, and the GM volume of the right inferior temporal gyrus, anterior division and the right accumbens.
For better illustration and clearer visualization of the difference for these selected features between PD group and HC group, in Table 3 we reported the mean, standard deviation (SD) and P-value of the resulting features from the dimension reduction step for two groups. We also plotted the histograms of FIGURE 1 | Flowchart of the study. After the rs-fMRI and sMRI images were preprocessed, we extracted the 6669 metrics. Then, Mann-Whitney U test, screening out high correlated variables and Lasso regression were implemented to reduce the number of features. Last, random forest and SVM were conducted to differentiate between PD and HC subjects.
these features with different colors representing different groups in Figure 3. In particular, the increasing or decreasing trends of these features between two groups can be immediately revealed by looking at the corresponding values in Table 3 or locations in Figure 3.

Model Fitting
After the screening process, we were left with only 56 features and were no longer stuck in the ultrahigh dimensional situation. Most of the existing machine learning methods including random forest and SVM could accommodate this relatively smaller number of variables compared with previous 6669 features. Therefore, we conducted the model fitting procedure for classifying subjects in the training set utilizing random forest and SVM to evaluate the performance of these two methods. It turned out that after cross-validation, both random forest and SVM achieved the perfect accuracy and AUC for distinguishing between PD and HC subjects in the training set, which was not surprising since we only had 56 coefficients to estimate while we had a sample size of 93 subjects.

Model Validation
Despite the superior performance of both methods in the training set, what really matters is the predictive result in the testing set. We therefore examined the validity of random forest and SVM by evaluating their classification performance in the testing set using AUC and the following measures (Table 4) TPR measures the proportion of PD patients correctly detected by the given procedure among all the PD patients.
FPR is calculated using the number of people who were falsely identified as having PD, divided by the total number of HCs. Accuracy gives the proportion of true results (both true positives and true negatives) among the total number of cases examined, i.e., the sample size of the testing set. In practice, we used the natural cutoff 0.5 to determine whether the subject should be classified as PD. To further assess the robustness of two methods, we also plotted the ROC curves in Figure 4 by varying thresholding values.
From Table 4 and Figure 4 we can tell that random forest and SVM performed comparably well in terms of the overall accuracy. Random forest performed slightly better than SVM according  The larger corresponding value on the color bar means more connections are involved in the respective region. The brain networks were visualized with the BrainNet Viewer (Xia et al., 2013).
to FPR, while SVM outperformed random forest according to accuracy, TPR and AUC. Therefore, for a less conservative and more robust prediction, we would prefer SVM over random forest based on the performance summary in the testing set.

DISCUSSION
We presented a framework for uncovering predictive markers of PD based on radiomics analysis and obtained excellent accuracy in classifying PD from HCs. Our model was established based on the relevant clinical characteristics, wholebrain functional connectivity and activity along with gray matter structure. After collecting the MRI scans for one subject, it may take 2 h for obtaining clinical evaluation, 10 min for preprocessing and around 2 s for running the machine learning models. To the best of our knowledge, this is the first study to explore the whole-brain functional activity and gray matter structure in a homogeneous and relatively large sample MRI study. The distinctive wholebrain functional activity and connection were mainly located within or across the AN, DMN, ECN and SMN in PD compared with HCs.
Our results showed that both methods achieved perfect accuracy in the training set, and SVM yielded an overall better classification performance than random forest in the testing set. In particular, SVM had higher accuracy (85%), TPR (1) and AUC (0.97) than random forest, while the FPR for SVM (0.31) was higher than random forest (0.27). The radiomics-based machine learning models in present study demonstrated the validity of trained classifiers in PD, which could be helpful to support clinical decision in both radiology and neurology.
Previous studies have made great progress in identifying PD and other neurodegenerative disorders from HCs using structural and fMRI data with the assist of machine learning. These results demonstrated the ability of supervised classification methods with a relatively high accuracy. An automatic SVM based study with leap motion controller recruited 16 PD and 12 HCs, and the accuracy was 74.07% with an AUC of 0.675 (Butt et al., 2018). Another study of SPECT imaging using SVM and logistic regression (LR) showed that SVM method produced a higher accuracy of 85% than LR of 83%, and the authors claimed that the SVM-based model could provide better guide for PD stage classification (Hsu et al., 2019). A large sample based on 831 structure T1-weighted MRI achieved a very high accuracy of up to 99% for differential diagnosis of 3 | The mean, standard deviation (SD) and P-value for all 54 selected features in PD group and HC groups.

ID
Feature PD (mean ± SD) HC (mean ± SD) P-value  PD (Singh and Samavedham, 2015). A study incorporating DTI and VBM in an SVM algorithm correctly distinguished PD from progressive supranuclear palsy (PSP) when white matter atrophy was considered (Cherubini et al., 2014). Combined with these previous findings, big data-driven approaches were helpful to aid PD diagnosis and to reach precision medicine (Khoury et al., 2019;van den Heuvel et al., 2020). The aforementioned methods either possessed a lower accuracy and AUC or only considered part of the complete radiomic features based on either rs-MRI or sMRI. In our study, the radiomics approach integrated both rs-MRI and sMRI by extracting features that quantify the whole-brain functional activity and connectivity along with GM volume and clinical evaluations. Forty-six RSFCs and the respective connected brain regions were selected after dimension reduction, and these disturbed brain regions related to RSFCs were primarily located in the ECN, DMN, AN, VN, and SMN. Seven more features of the HAMD, mALFF, mReHo, VMHC and the GM volume were also retained to build the classifying model. Within the model, intrinsic connectivity networks of mALFF were identified and located mainly in DMN including the left superior temporal gyrus, posterior division, the left parahippocampal gyrus, posterior division. The particular neural activities of mReHo maps were also located in DMN of the right thalamus and GN of the left pallidum. Selected VMHC quantified the functional homotopy of connectivity in VN of the right temporal fusiform cortex and anterior division. Features of the GM volume also covered the DMN in the right inferior temporal gyrus, anterior division and the right accumbens. Some of these findings were expected and in accordance with previous studies. For example, the altered RSFCs were primarily located in the typical resting-state network (RSN). However, the prominent role of DMN, ECN, AN, VIN and sensorimotor functioning in PD revealed in our study was remarkable. The identified regions in DMN included the left superior temporal gyrus, posterior division, the left parahippocampal gyrus, the right thalamus and the right inferior temporal gyrus. We also found abnormal AN in the right insular cortex, left and right temporal pole, along with unusual VIN in lateral occipital cortex, superior, temporal fusiform cortex, the left and right temporal occipital fusiform cortex, the left occipital pole, the right superior parietal lobule. Abnormal ECN was detected in the left and right paracingulate gyrus, and the right anterior cingulate gyrus.
RSN reflects the spontaneous neural activities of the blood oxygenation level-dependent (BOLD) signals between temporally correlated brain regions. Compared with the control group, the DMN plays a crucial role in neurodegenerative disorders and normal aging. Several fMRI studies have indicated that the DMN injured before the cognitive decline in PD (Sandrone and Catani, 2013;Koshimori et al., 2016). A 2-year study using ReHo and VBM to identify differences in local spontaneous brain activity and gray matter volume found that PD patients with normal cognition showed a decreased ReHo in the DMN (Zeng et al., 2017). In addition, a gender-specific effect of uric acid on restingstate cortical FC found the de novo PD group had decreased FC in bilateral cingulate, postcentral and lingual gyri within DMN (Lee et al., 2018). Our results were consistent with these previous studies.
The basal ganglia, thalamus, and brainstem are important in the pathophysiology of PD. Studies on detecting neural activity changes in these regions have achieved more sensitive and reliable results for scientific and clinical research on PD. The basal ganglia network (BGN) has been observed in pathologies with altered neurotransmitter systems of dopaminergic processes, and also involving motor control. In the present study, we found disturbed BGN in left pallidum, the right thalamus and the right brainstem. The variability of FC in healthy older adults found strongest correlate of FC in the BGN, and potential links to dopaminerelated function (Griffanti et al., 2018). A sex-related pattern of RSN showed an increased connectivity within the BGN in female PD patients, and FC changes in sensorimotor at baseline were considered as an independent predictor of disease severity in the early stage of PD (De Micco et al., 2019).
Attention networks (AN) in cortical regions are affected in early stage of PD (Madhyastha et al., 2015). Proteinopathy and longitudinal changes in FC networks within the SMN were confirmed, and the interaction between the dorsal attention network (DAN) the frontoparietal control network decreased significantly over time in PD while correlated with the decline in cognitive function (Campbell et al., 2019). Altered organization of the DAN and lack of changes in the ventral attention network (VAN) in PD patients indicated the higher risk for freezing of gait during complex walking situations, and these findings revealed that AN played an important role in freezing of gait (Maidan et al., 2019). Gender-specific effect of uric acid on rs-fMRI networks in de novo PD found decreased FC in bilateral insular, frontal and temporal areas within DAN and bilateral medial temporal and right insular areas within executive control network (ECN) (Lee et al., 2018).
Apart from these RSFC findings, structure MRI has received more research focus on better stability and repeatability compared with rs-fMRI. In our study, GM volume of the right inferior temporal gyrus, anterior division and the right accumbens demonstrated differences in PD as opposed to HC. Atrophy of the putamen and altered FC of the striatal structures in PD revealed key structure-function relationship. Caudate nucleus and putamen atrophy could serve as neuroimaging biomeasures for PD (Owens-Walton et al., 2018). A brain microstructual study found decreased white matter fiber features in the right arcuate fasciculus and bilateral middle cerebellar peduncles. The study also detected increased network connectivity in prodromal early PD, which might indicate the neural compensation (Sanjari et al., 2019). The right accumben as one of selected GM features in the model is an interesting sign. An analysis of dopamine regulation and transporter function found regional brain (the nuclei accumbens, cingulate regions and inferior frontal) were closely related with apathy rating scores and β-amyloidopathy for predicting cognitive decline in advancing PD (Zhou et al., 2020). An eventrelated fMRI study based on reward-related neural responses showed the left nucleus accumbens with lower activation indicated involvement of the ventral striatum in individuals for further development of PD (Thaler et al., 2019) and PD patients with persistent pain displayed an accumbens-hippocampus disconnection (Polli et al., 2016).
Our study certainly has several limitations. First, although we have carried out detailed clinical evaluation and stage classification for all subjects, due to our limited sample size, we could not further stratify the patients according to the disease severity. Second, RSFCs in the model may be influenced by the different clinical symptoms such as olfaction or depression. Third, cerebellum networks were left out in the model. Although the key molecular events that provoke PD have not been fully understood and the underlying mechanisms involving cerebellum were relatively less reported, considerable evidence has indicated that cerebellum plays an important role in sensorimotor dysregulation and has now received growing attention (Lopez et al., 2020;Maas et al., 2020). For future studies, we will include the cerebellum, and increase our sample size in order to obtain different subgroups based on disease stages.
In agreement with previous rs-fMRI studies, the proposed radiomics method that combined rs-fMRI spontaneous activity, connectivity and structure MRI of gray matter (GM) was proved to be scientifically sound and valid. The machine learning based radiomics approach can help the diagnosis, personalized treatment, and prognosis orientation for patient with PD at a lower cost. This type of radiomics approaches should be widely performed and considered as an automated classification framework for predicting PD in clinical management.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Medical Research Ethical Committee of Nanjing Brain Hospital (Nanjing, China). The patients/participants provided their written informed consent to participate in this study.