Convolutional Neural Network Optimize the Application of Diffusion Kurtosis Imaging in Parkinson’s Disease

Junyan Sun Department of Neurobiology, Neurology and Geriatrics, Xuanwu Hospital of Capital Medical University, Beijing Institute of Geriatrics,Beijing, China https://orcid.org/0000-0002-2524-6133 Ruike Chen Center for Brain Imaging Science and Technology, College of Biomedical Engineering and Instrumental Science, Zhejiang University, Hangzhou, Zhejiang,China. 310027 Qiqi Tong Research Center for Healthcare Data Science, Zhejiang Lab, Hangzhou, Zhejiang,China Jinghong Ma Department of Neurology, Xuanwu Hospital of Capital Medical University,Beijing,China Linlin Gao Deparment of Neurology, Neurology and Geriatics, Xuanwu Hospital of Capital Medical University, Beijing Institute of Geriatrics, Beijing, China Jiliang Fang Department of Radiology, Guang'anmen Hospital, China Academy of Chinese Medical Sciences, Beijing,China Dongling Zhang Department of Neurology, Neurology and Geriatrics, Xuanwu Hospital of Capital Medical University, Beijing Institute of Geriatrics, Beijing, China Piu Chan Department of Neurobiology, Neurology and Geriatrics, Xuanwu Hospital of Capital Medical University, Beijing Institute of Geriatrics, Beijing, China Hongjian He Center for Brain Imaging Science and Technology, College of Biomedical Engineering and Instrumental Science, Zhejiang University, Hangzhou, Zhejiang, China. 310027 Tao Wu (  wutao69@gmail.com ) Department of Neurobiology, Neurology and Geriatrics Xuanwu Hospital of Capital Medical University


Introduction
Parkinson's disease (PD) is a common neurodegenerative disease characterized by bradykinesia, resting tremor, rigidity, posture balance disturbance and non-motor manifestations 1 . Beyond the de ciency of dopaminergic neurons and aggregation of Lewy bodies in the basal ganglia, pathological changes of PD are associated with axonal lesions and synaptic dysfunction, which contributes to the impairments of white matter integrity 2 . Given the limitation of discerning the intrinsic details and pathological heterogeneity in brain tissues, conventional magnetic resonance imaging (MRI) has di culty in revealing PD-associated microstructural changes. Diffusion-weighted MRI techniques like diffusion tensor imaging (DTI), can non-invasively probe the microstructural properties via the diffusion of water molecules in vivo [3][4][5][6][7] . It has been reported that DTI-derived metrics, such as fractional anisotropy (FA) and mean diffusivity (MD), showed signi cant differences in the substantia nigra (SN) and some white matter areas in PD patients 8 . Further, diffusion kurtosis imaging (DKI), which developed based on DTI and taken the non-Gaussian diffusion of water molecules into account, was reported as a more sensitive technique to evaluate the pathological characteristics of PD patients 5,8 .
However, previous studies yielded inconsistent or even controversial ndings. For example, while some studies showed decreased FA value, increased MD and/or mean kurtosis (MK) values in the SN 6,[9][10][11] , other studies 6, 12 observed increased FA value in the SN in PD patients. It has been reported that the FA value of SN in PD patients had a rising tendency compared with healthy controls (HCs) 11 . Additionally, Kamagata  It is speculated that the heterogeneity of PD patients being recruited and various acquisition protocols of diffusion MRI scanning may contribute to those controversial ndings 3,4,8,14 . DKI can sensitively re ect microstructural complexity, particularly in isotropic tissues such as gray matter 15 . But since gray matter microstructure lacks evident directionality, DWI signals can be easily affected by noise and limited spatial resolution 16 , thus leading to inaccurate ndings of alterations in DKI scalar measures.
Recently, deep learning, an important branch of machine learning, shows noteworthy potential for improving the detection of neuroimaging ndings [17][18][19][20][21] . As one of the representative algorithms of deep learning, convolutional neural network (CNN) adopts convolution and down sampling to certain layers with less computation, adjusts the network weights through the back-propagation and stochastic gradient descent algorithm, recognizes the features or patterns of the raw imaging inputs automatically, and then achieves the classi cation, identi cation and prediction of inputs 19,21,22 .
Li et al. 23 recently proposed a three-dimensional hierarchical convolutional neural network (3D H-CNN) to improve the estimation of DKI scalar measures from limited diffusion-weighted images (DWIs). Threedimensional convolution kernels were introduced to extract and learn features of the DWIs automatically. Only part scalar measures were of clinical interest instead of the full tensors, 3D H-CNN (herein after called CNN) method makes it possible to complete fast and optimize DKI acquisition within 1 minute. It also took the cross-voxel information into account, which was proved to provide enhanced e ciency of estimating DKI scalar measures and improved robustness against noise.
Therefore, in the current study, we aimed to use this CNN method to improve the estimation of DKI scalar measures, and to examine whether the improved measures can help reveal PD-associated imaging features.

DATASET 2
Five travelling subjects were scanned on two different MRI scanners. Scan 1 was also using the scanner-A and same scanning process as performed in the DATASET 1. For implementation of CNN, DWIs were also acquired in the Scan 2 using a 3T scanner-B (MAGNETOM Prisma, Siemens Germany) equipped with a 64-channel RF coil. The DWIs were obtained using a simultaneous multi-slice diffusion echo planar imaging sequence (SMS-EPI). Diffusion weightings of b = 1000, 2000 and 3000 s/mm² were applied in 30 gradient directions. Six b = 0 images and one b = 0 image with an opposite phase encoding direction were acquired,resulting in a total of 96 DWIs. The one b = 0 image with reversed phase encoding direction was for correcting the eld inhomogeneity-induced distortion. Other imaging parameters were as follow: repetition time (TR) = 3000 ms, echo time (TE) = 75 ms, resolution = 2mm×2mm×2mm, FOV = 220mm×220 mm, number of slices = 68.

preprocessing
The preprocessing pipeline for both datasets was mainly based on FSL (FMRIB Software Library, University of Oxford, UK) 24 . Rician noise was removed using dwidenoise included in MRtrix 3 25 , followed by gibbs-ring removal. For DATASET 2, Scan 2, distortion correction was also performed using topup and eddy tool 26 . The susceptibility-induced off-resonance eld map was rst estimated by topup using a pair of non-DW (b = 0) images acquired with reversed phase encoding directions AP and PA. It was then fed into eddy to perform correction for eddy current and motion induced distortion.

model t
For the DATASET 1 and Scan 2 of DATASET 2, model tting method was conducted using DESIGNER (Diffusion parameter EStImation with Gibbs and NoisE Removal, New York University, US), which is a post-processing pipeline capable of identifying and correcting various speci c artifacts and confounding factors for an improved accuracy, precision, and robustness compared to traditional linear least square method tting ( Figure. 1) 27 .

CNN
CNN method was adopted to improve the estimation of DKI scalar measures from limited quality of DATASET 1. The adopted CNN method included one input layer, several hidden layers and two output layers. A dropout layer was inserted before each output layer to prevent over tting. 3×3×3 convolution kernels were introduced in the rst hidden layer to extract features from the preprocessed DWIs. The network was constructed in a hierarchical structure. Resulting DKI scalar measures were output through two different layers. The shallow output layer was connected to the penultimate hidden layer and was responsible for the scalar measures (FA and MD). Kurtosis related measures (MK and KFA) values were output through a deeper layer connected to the last hidden layer 23 .
Preprocessed DWIs of DATASET 2, Scan 1 were the training data of CNN and the corresponding modeltted MK, KFA, FA, and MD maps for each traveling subject in DATASET 2, Scan 2 were de ned as training and test labels. Scan 2 were registered to Scan 1 using non-linear registration before training. Preprocessed DWIs of DATASET 1 were used as test sets ( Figure. 1).
The entire training and testing pipeline of CNN is as displayed in Fig. 2. Preprocessed DWIs of Scan 1 and the reference standard DKI maps derived from scan 2 were input as the training and labeling datasets in the training phase. Noise and gibbs-ring removed DWIs of DATASET 1 were the inputs as the testing phase. DWIs of Scan 2 included higher b values and more diffusion weighted directions than Scan 1. The additional b = 0 images with reversed phase encoding directions allows correction for eld inhomogeneity induced distortion. All improvements in scanning condition resulted in a better quality of DWIs and a better estimation of DKI scalar measures compared with model tting method. The CNN training process related limited quality DWIs with its corresponding high-accuracy DKI scalar measures, gaining the network an ability to predict high-quality DKIs from moderate DWIs ( Figure. 1). As shown in Figure. 1 and Figure. 2, DKI scalar maps estimated by CNN of both dataset 1 and 2 had higher signal-to-noise ratio (SNR) than those by model tting. For dataset 2 the CNN results displayed almost equivalent quality as the reference standard (ground truth).

Statistical analysis
The 1-sample Shapiro-Wilk test was applied to con rm the normality of each group data. The Student's t test and Pearson chi-square test were used for age and sex variables, respectively.

Whole-brain unpaired t-test
For the DWIs of PD patients and HCs, two datasets of DKI scalar measures were derived using the model tting method and CNN method, respectively. FA maps were rst registered to the MNI 152 standard space with resolution of 1mm isotropic using a combination of linear and non-linear transforms. The resulting transformation was then applied to all other DKI maps for co-registration. For DKI scalar measures of both methods, whole-brain unpaired t-tests were performed between PD patients and HCs to evaluate the ability of DKI to reveal the differences between the groups.

Tract-Based Spatial Statistics
Furthermore, we performed the Tract-Based Spatial Statistics (TBSS) analysis using the model tting method and CNN method, respectively. TBSS extracted the mean FA maps to generate white matter skeleton, realized by a tool for nonparametric permutation inference implemented in FSL 28 .
Threshold-Free Cluster Enhancement (TFCE) 29 based test was included in both whole-brain unpaired ttest and TBSS analysis to improve robustness compared with conventional voxel-based tests. All signi cance threshold was determined at p < 0.05 and corrected by family wise error (FWE).
It is noteworthy that we validated whether there were intra-group differences in HCs applying the analysis of whole-brain unpaired t-test and TBSS with model tting method and CNN method respectively before comparing differences between the groups.

Correlation analysis
To determine the clinical signi cance of DKI scalar measures with CNN method more clearly, DKI scalar measures (MK, KFA, FA and MD) with model tting or CNN method showing signi cant between-group differences in the basal ganglia (SN, putamen and caudate nucleus) were extracted and correlated with clinical assessments. Pearson's correlation analysis was used for normal distribution data and Spearman's correlation was used for non-normal distribution data. Correlation with signi cance was de ned at p < 0.016 (Bonferroni corrected). Statistical analyses were computed using IBM SPSS Statistics (version 25; IBM Corp., Armonk, NY) and GraphPad Prism 8.0.1.

Demographic features
Age (p = 0.659, Student's t test) and sex distribution (p = 0.092, Pearson chi-square test) did not differ between PD and controls (Table 1).

The validation of intra-group difference in HCs
There were no signi cant differences of DKI scalar measures within healthy controls with model tting and CNN methods.

model tting method
FA values in the bilateral anterior and superior corona radiata increased in PD patients compared to HCs. MD values in the caudate, putamen and thalamus bilaterally, and bilateral cerebral cortex and white matter, bilateral posterior thalamic radiation (include optic radiation), right sagittal stratum (including inferior longitudinal fasciculus and inferior fronto-occipital fasciculus), right cingulum (hippocampus), right superior and anterior corona radiata, and left fornix were signi cantly increased in PD patients compared with HCs. There was no signi cant difference of MK and KFA values between the groups ( Figure. 3) (p < 0.05, FWE corrected).

Discussion
The major nding of this study was that the CNN-estimated MD values in the left SN and bilateral caudate increased in PD patients compared to HCs. Additionally, the CNN-estimated FA and MK values in the right SN negatively correlated with H&Y scales and CNN-estimated FA values in the left putamen positively correlated with H&Y scales. In contrast, with model tting method, there was no signi cant difference of MD values in the SN between PD patients and HCs, and there was no signi cant correlation between DKI scalar measures and clinical assessments in PD patients. Our ndings suggest that the CNN method have the potential to optimize the estimation of DKI scalar measures and to improve the sensitivity to detect PD-related imaging features.
With CNN method, we found greater MD values in several brain regions in PD patients compared with HCs, especially in the left SN, which is consistent with previous reports using regions of interest (ROI) analysis 10,11,14,[30][31][32] . PD is characterized by progressive death of dopaminergic neurons in the SN, followed by the loss of dopaminergic projections from the SN to striatum, resulting in a series of motor and non-motor symptoms 33 . According to the mathematical concept of tensor, the three-dimensional shape of diffusion elliptical structure depends on three eigenvalues (λ1, λ2, λ3) of orthogonal principal axes without directions. The MD value is the average of the three eigenvalues. The impaired axons and neurons and loss of myelin integrity in PD patients result in the decrease of restriction of water molecules displacement, which induces increased MD values 32, 34 . Regional increased MD values in the left SN and bilateral caudate estimated by CNN method are consistent with the pathological lesions in PD patients. In contrast, we did not nd increased MD in the SN in PD patients compared to HCs by applying model tting method. This nding indicates that CNN method can better reveal the pathological features of PD than model tting method.
We did not observe modulation of FA and MK values in the SN in PD patients, which is in line with a previous report 14 . In contrast, some previous studies based on ROI analysis showed decreased or increased FA and/or increased MK in the SN in PD patients 6,31 . We speculate that, rst, different analysis methods may responsible for those controversial results. Whole-brain unpaired t-test, moving beyond the hypothesis-driven ROI analysis, focused the statistical information on each voxel accompanied with increased partial volume effects and false-positive risk, especially within the pathological brain tissues. Second, we suppose that these controversial ndings may due to the heterogeneity of patients being recruited and variations imaging quality 3,8 . In addition, it has been reported that iron deposition could increase FA values and decrease MD values in the white and gray matters 35 . Numerous reports have demonstrated that there was iron accumulation in the SN [36][37][38] . Thus, different levels of iron deposition in the SN may also be a reason contributing to these inconsistent ndings.
We found negative correlation between H&Y scales and CNN-estimated FA values and MK values in the SN, as well as positive correlation between H&Y scales and CNN-estimated FA values in the putamen.
These results indicate that FA and MK in the SN decreases, while FA in the putamen increases with the progression of the disease. As most of our patients were in the early stages ( fty-ve of our patients were at H&Y stages 1 and 2), it is possible we could detect decreased FA in the SN if more advanced patients were enrolled. We did not nd any signi cant correlation between DKI scalar measures and clinical assessments in PD patients using the model tting method, which further proves that using CNN method to estimate DKI measures can improve the ability to explore PD-related neural modulations than using model tting method.
For the TBSS analysis, increased FA values were observed in the brain white matter such as anterior thalamic radiation, inferior longitudinal fasciculus, superior longitudinal fasciculus, corticospinal tract, and inferior fronto-occipital fasciculus with both methods and agreed with previous studies 13, [39][40][41] . It has been shown that increased FA in these white matters correlated with better olfactory function and lower motor severity 42 . Thus, these increased diffusional properties of white matter might be a re ection of microstructural compensation 42 .
We observed greater MK values in white matter in PD patients, which is inconsistent with previous reports. Previous studies found either no signi cant difference in MK values 43,44 , or decreased MK values of anterior cingulum, inferior fronto-occipital fasciculus and uncinate fasciculus in PD patients 7,45 . We suppose that the heterogeneity of patients being recruited and differences in protocol of diffusion image acquisition and image processing may contribute to these inconsistent ndings. In addition, we found increased KFA in white matter, which has not been reported previously. KFA values, resembling the FA de nition, quantify the degree of anisotropy of the non-Gaussian diffusion. In the current study, the increased KFA and FA values are located in the same white matter bers. So far, only a small number of studies have paid attention to the kurtosis changes of white matter in PD patients, it is necessary to perform large cohort studies to elucidate the microstructural changes in white matter in PD patients.
In conclusion, the CNN method has the potential to sensitively detect the nigral pathology and improve the robustness and performance of DKI images with few of DWIs and then to differentiate PD patents from HCs.

Declarations
Availability of data and materials The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Competing interests
Junyan Sun and Ruike Chen analyzed the data and wrote the manuscript. Qiqi Tong guided the analyses of diffusional MRI data. Tao Wu and Hongjian He designed and supervised the whole research and polished the manuscript. Jinghong Ma, Linlin Gao, and Dongling zhang collected clinical assessment of the recruited subjects. Jiliang Fang is responsible for the scanning and quality control of raw MRI data.
Piu Chan polished the manuscript. All authors read and approved the nal manuscript." Figure 1 The training and testing pipeline of CNN method. Note. CNN=convolutional neural network; DWI=diffusion weighted image; DKI=diffusion kurtosis imaging; FA=fractional anisotropy; MD=mean diffusivity; MK=mean kurtosis; KFA=kurtosis fractional anisotropy.

Figure 3
The whole-brain unpaired t-test analysis of DKI measures with CNN and model tting methods. Note.
Increased MD values in bilateral caudate (blue circle) and right SN (yellow circle) with the CNN method between HC and PD groups (p< 0.05, corrected by family wise error). KFA=kurtosis fractional anisotropy; FA=fractional anisotropy; MD=mean diffusivity; CNN= convolutional neural network; HC=healthy controls; PD, Parkinson's disease; MNI=Montreal Neurological Institute coordinate.