Pancreatic Mass Characterization Using IVIM-DKI MRI and Machine Learning-Based Multi-Parametric Texture Analysis

Non-invasive characterization of pancreatic masses aids in the management of pancreatic lesions. Intravoxel incoherent motion-diffusion kurtosis imaging (IVIM-DKI) and machine learning-based texture analysis was used to differentiate pancreatic masses such as pancreatic ductal adenocarcinoma (PDAC), pancreatic neuroendocrine tumor (pNET), solid pseudopapillary epithelial neoplasm (SPEN), and mass-forming chronic pancreatitis (MFCP). A total of forty-eight biopsy-proven patients with pancreatic masses were recruited and classified into pNET (n = 13), MFCP (n = 6), SPEN (n = 4), and PDAC (n = 25) groups. All patients were scanned for IVIM-DKI sequences acquired with 14 b-values (0 to 2500 s/mm2) on a 1.5T MRI. An IVIM-DKI model with a 3D total variation (TV) penalty function was implemented to estimate the precise IVIM-DKI parametric maps. Texture analysis (TA) of the apparent diffusion coefficient (ADC) and IVIM-DKI parametric map was performed and reduced using the chi-square test. These features were fed to an artificial neural network (ANN) for characterization of pancreatic mass subtypes and validated by 5-fold cross-validation. Receiver operator characteristics (ROC) analyses were used to compute the area under curve (AUC). Perfusion fraction (f) was significantly higher (p < 0.05) in pNET than PDAC. The f showed better diagnostic performance for PDAC vs. MFCP with AUC:0.77. Both pseudo-diffusion coefficient (D*) and f for PDAC vs. pNET showed an AUC of 0.73. ADC and diffusion coefficient (D) showed good diagnostic performance for pNET vs. MFCP with AUC: 0.79 and 0.76, respectively. In the TA of PDAC vs. non-PDAC, f and combined IVIM-DKI parameters showed high accuracy ≥ 84.3% and AUC ≥ 0.84. Mean f and combined IVIM-DKI parameters estimated that the IVIM-DKI model with TV texture features has the potential to be helpful in characterizing pancreatic masses.


Introduction
Pancreatic cancer is a challenging malignancy with a 5-year survival rate of only 10% due to its non-specific symptoms at an early stage [1]. Pancreatic masses are a spectrum of tumors, with pancreatic ductal adenocarcinoma (PDAC) accounting for 90% of all cases and having the worst prognosis. In contrast, pancreatic neuroendocrine tumors (pNETs), mass-forming chronic pancreatitis (MFCP), cystic masses, and solid papillary epithelial neoplasms (SPENs) account for less than 10% of pancreatic tumors with relatively good prognoses [2][3][4]. On imaging, PDAC, MFCP, and pNET may appear similar and distinguishing between these masses with currently available imaging modalities is difficult.
Several diagnostic tools, such as the serum levels of carbohydrate antigen (CA , have been used to characterize pancreatic masses [5]. However, they have a high falsenegative rate, and at the same time benign pancreatic masses can also have high CA 19-9 levels [5]. Fine needle aspiration cytology (FNAC) produces reliable results, but it is susceptible to tumor seeding in the sampling tract and has limitations such as sampling errors [6]. Cystic masses are easily distinguished by their morphological appearance on cross-sectional imaging with CT and MRI [7], however there are phenotypic similarities between other benign entities like MFCP and malignancies like PDAC [8].
Non-invasive assessment of solid pancreatic masses is preferable, and it is feasible with diffusion-weighted imaging (DWI), which is widely recognized for characterization of tumor physiology using apparent diffusion coefficients (ADCs) [9]. The ADC, which is derived from the DWI signal, has been found to capture the change in tumor cellularity [9]. However, ADC is sensitive towards interfering perfusion signals, which becomes significant as the diffusion gradient strength decreases (at lower b-values) [10,11]. Intravoxel incoherent motion (IVIM) can be useful to study the perfusion information separately from the acquired diffusion data [10]. However, this model assumes a Gaussian distribution of underlying water movement, which may result in estimation errors for diffusion parameters [12]. This is not the case with diffusion kurtosis imaging (DKI), which corrects for the non-Gaussian distribution of water diffusion caused by underlying tissue microstructures [12]. Lu et al. proposed simultaneous modeling of diffusion (molecular diffusion (D)) and perfusion parameters (pseudo-diffusion coefficient (D*) and perfusion fraction (f)) with a kurtosis (k) measure of non-Gaussian distributions of water movements using a hybrid IVIM-DKI model [13]. This model accurately captures the tumor heterogeneity as compared to other DWI models, allowing improved tumor characterization as shown in brain [14], head and neck [13], thyroid [15], and prostate [16] tumors.
The existing research gap in the existing methodology is as follows: simultaneous modelling of all parameters might result in non-physiological heterogeneity in the parametric maps, possibly explaining the lack of wide acceptance of this imaging modality in routine clinical use for a long time. Table 1 shows the advantages and drawbacks of existing techniques. Total variation penalty function (TV) can be used with non-linear least-square optimization; this incorporates a physiologically reasonable spatial constraint into the standard IVIM-DKI model. Previous studies have demonstrated that this method was useful for the removal of non-physiological heterogeneity in the parametric map by iteratively lowering spurious values throughout the parameter estimation [17,18].

DWI/ADC
DWI allows for the qualitative and quantitative evaluation of tissue diffusivity without the need for contrast agents. ADC is quantified using a simplistic model and provides quantifiable measures of tumor cellularity.
ADC can be affected by perfusion signals from flowing blood effects that can cause overestimations. ADC becomes more sensitive to tissue microscopic characteristics with higher b-values due to non-Gaussian water diffusion.

IVIM/D/D*/f
Diffusion and perfusion components of tissue can be evaluated independently.
D* and f can be overestimated due to non-Gaussian water diffusion at high b-values and long scan times.
DKI/D/k Accurate representation of water interactions inside tumor.
DKI model cannot eliminate the perfusion effect at low b-values.
Low SNR with noisy parametric maps and long scan times.
Despite developments in CT and MRI technologies, detecting and characterizing pancreatic masses remain difficult. Texture analysis (TA) has shown promise in early diagnosis by collecting imaging characteristics that offer connections between quantitative biomarkers, therapeutic response, and surgical planning [19]. By calculating gray-level patterns and local spatial correlations between patterns and textures in the area, the extracted features give information on inherent characteristics of tumor physiology and heterogeneity [20,21]. Only a few studies have utilized TA on MRI in the characterization of pancreatic masses; it has considerably assisted in the distinction of pancreatic masses [22][23][24][25]. Thus, the main objective of this study was to evaluate the added value of IVIM-DKI analysis with a TV penalty function in texture analysis in the characterization of pancreatic lesions. Second, we evaluated the role of multi-parametric texture analysis of IVIM-DKI for the characterization of pancreatic masses utilizing machine learning-based techniques. The main contributions of this study are:

•
We established a novel IVIM-DKI model with a total variation penalty function to achieve improved non-invasive characterization of pancreatic masses.

•
Qualitative and mean comparison between IVIM-DKI parametric maps in pancreatic masses such as PDAC, PNET, MFCP, and SPEN were evaluated. • Cut-off values for each IVIM-DKI parameter were calculated for characterization of pancreatic masses using ROC analysis.

•
We attempted to comprehensively investigate texture features of apparent diffusion coefficient (ADC), diffusion coefficient (D), pseudo-diffusion coefficient (D*), perfusion fraction (f), kurtosis (k), and combined texture features of IVIM-DKI parameters with and without ADC. • Machine learning-based classification of pancreatic masses using ANN was used and compared with other techniques such as decision tree and ensemble.

Study Population
Patients were recruited for this prospective study after approval from the institute ethics committee and informed written consent was obtained from the study participants before enrollment. Patients were recruited and categorized as having a solid or cystic pancreatic mass based on preliminary MRI. A total of seventy-eight patients (n = 78) with pancreatic masses were included from May 2017 to June 2019. Enrolment exclusion criteria of this study were patients with: 1. unequivocal pseudocyst; 2. collection in acute/chronic pancreatitis; 3. prior history of treatment/surgery for a pancreatic mass; 4. contraindications to MRI; 5. terminally ill-patients; and 6. refusal of consent. However, only forty-eight patients could be included in this study analysis out of total of seven-eight patients because the imaging workup was incomplete in six patients (n = 6), the mass was of extra-pancreatic origin on histopathology analysis in four patients (n = 4), the mass was a complication of acute pancreatitis in ten patients (n = 10), and ten patients had unequivocal cystic masses (n = 10). A final sample size of forty-eight (n = 48) was used for further analysis. Pancreatic masses were proven either on histopathology of the surgically resected specimen or on fine needle aspiration cytology (FNAC) of the mass and they were considered as reference standards for thirty-eight (n = 38) patients and for remaining ten patients (n = 10) imaging was considered as the gold standard. FNAC was performed by using either endoscopic ultrasound (EUS) guidance or transabdominal ultrasound guidance.

MRI Acquisition
All patients underwent MRI on a 1.5T scanner (Achieva, Philips Healthcare, Best, The Netherlands) using a multi-channel phased-array body coil. The MRI sequence protocol used to acquire all the patients' scans included gradient echo T1-weighted (mDixon) with repetition time (TR): 500 ms and echo time (TE): 2.3 and 4.6 ms, fat-suppressed (FS) T2-weighted (turbo spin-echo (TSE)) was acquired in the axial plane with TR: 1000 ms and TE: 80 ms. Steady-state free precession (SSFP) sequences were acquired in the coro-nal plane with TR: 500 ms and TE: 50 ms. Thick slab half-Fourier acquisition single-shot turbo spin echo (HASTE) with TR: 8000 ms and TE: 800 ms and driven equilibrium with  90-degree flip-back pulse (RESTORE) with TR: 1000 ms and TE: 6500 ms sequences for  heavily T2-weighted Thick slab magnetic resonance cholangiopancreatography (MRCP) sequence (MRCP) were performed. Contrast-enhanced T1-weighted mDixon sequences were acquired in the arterial, pancreatic, and venous phases at 25 s, 45 s, and 70 s, respectively. The patient received 0.2 mL/kg of intravenous gadolinium-based MRI contrast by a pressure injector, followed by a flush of saline.
IVIM-DKI sequences were acquired using free-breathing technique with TR = 1000 ms, TE = 118.7 ms, slice thickness = 6 mm, the field of view (FOV) = 375 × 305 mm, and matrix size = 124 × 100. IVIM-DKI images were acquired at 14 b-values covering the pancreas in the transverse plane using spin echo-echo planar imaging (SE-EPI) sequence with b-values (number of averages) of 0 , and 2500(7) s/mm 2 . The diffusion sensitizing gradient was applied in three orthogonal directions. Three raw IVIM-DKI images were acquired (one for each diffusion direction) and reconstructed. Then, these three images were averaged in the console to generate one IVIM-DKI image. The total scan time to acquire the IVIM-DKI images was 16 min per patient.

MRI Image Analysis
All IVIM-DKI images were exported from the picture archiving and communication system (PACS) workstation to a personal workstation in digital imaging and communications in medicine (DICOM) format and then converted to neuroimaging informatics technology initiative (NIfTI) format using the dcm2nii tool [26]. Non-linear least-squares optimization and parallel computing with an in-house built toolbox were used for ADC estimation. The ADC parametric map was estimated voxelwise using a monoexponential model with 4 b-values (i.e., 0, 500, 800, 1000 s/mm 2 ) as used in the literature [27,28] and in Equation (1): where S 0 and S indicate signal intensities of the images without and with diffusion weighting of b in s/mm 2 . The IVIM-DKI signal was modelled using two models with 14 b-values (i.e., 0 to 2500 s/mm 2 ): (1) standard model (Equation (2)), where IVIM-DKI parametric maps were obtained by simultaneously fitting for all the four IVIM-DKI parameters (D, D*, f, and k); and (2) IDTV model, where the IVIM-DKI model with TV method was used to estimate IVIM-DKI parametric maps using the open access ivimDKI3Dtvtool_v1_4 toolbox developed in MATLAB (version 9.9, The MathWorks, Inc., Natick, MA, USA) [16]. A trust-region based fit algorithm was used to estimate the IVIM-DKI parameters which allows us to use the bounds for each parameter, D = [0 0.05] mm 2 /s, D* = [0 0.5] mm 2 /s, f = [0 1], and k = [0 3] and was initialized with D = 1.3 × 10 −3 mm 2 /s, D* = 13 × 10 −3 mm 2 /s, f = 0.3, and k = 0.7 based on previous studies [29,30].
The IVIM-DKI images at different b-values for one representative patient are presented in Supplementary Material Figure S1. All IVIM-DKI signals were normalized to the b-value = 0 s/mm 2 signal and Equation (2) shows an association between signal variation and b-values in the IVIM-DKI sequence [13]: The IDTV model can remove any abrupt change in parameter estimates during the optimization process. The TV penalty function was used to iteratively balance out the estimated parameter [18]. The TV of a 3D parametric map is a L 1 norm of discrete gradient of the map [16]. The three dimensional-TV (3D-TV) is an adaptive regularization approach that has the benefit of being robust to b-values combinations and consistently improved parameter estimations at different field strengths [16,31]. The 3D-TV was used in conjunc-tion with the IVIM-DKI model to reconstruct the entire parametric map by considering neighboring voxels in three directions i.e., TV minimization was applied in the entire image at once [16].

Localization of Region of Interest
A radiologist with more than three years of expertise in abdominal imaging performed the region of interest (ROI) analysis, which was confirmed by another radiologist with more than ten years of experience in abdominal MR imaging and blinded to the reference standard. A free hand ROI was drawn over the lesion to cover the whole mass using MRICron software [26], as shown in Figure 1. The signal intensity of the mass was noted on T1-weighted, T2-weighted, DWI, and ADC images. As shown in Figure 1, the tumor ROI was drawn on the solid margin of the pancreatic tumors as seen on IVIM-DKI at b = 0 s/mm 2 and ROI was verified by the presence of hyperintensity and hyperintensity on corresponding slices of IVIM-DKI at b = 2500 s/mm 2 and ADC map, respectively. Mass ROIs were defined on all matching slices of IVIM-DKI images, using ADC maps and conventional sequences such as T2-weighted, contrast enhanced T1-weighted, and DWI images as a reference. age at once [16].

Localization of Region of Interest
A radiologist with more than three years of expertise in abdomin formed the region of interest (ROI) analysis, which was confirmed by an with more than ten years of experience in abdominal MR imaging and bl erence standard. A free hand ROI was drawn over the lesion to cover using MRICron software [26], as shown in Figure 1. The signal intensity noted on T1-weighted, T2-weighted, DWI, and ADC images. As shown tumor ROI was drawn on the solid margin of the pancreatic tumors as se at b = 0 s/mm 2 and ROI was verified by the presence of hyperintensity an on corresponding slices of IVIM-DKI at b= 2500 s/mm 2 and ADC map, re ROIs were defined on all matching slices of IVIM-DKI images, using ADC ventional sequences such as T2-weighted, contrast enhanced T1-weighte ages as a reference.    Figure 2 represents the flowchart of the machine learning-based classification of pancreatic masses using texture features from ADC and IVIM-DKI parameters. Three-dimensional TA was performed for the classification of pancreatic masses using the whole volume ROI which was superimposed onto ADC and IVIM-DKI parametric maps. A total of 30 texture features were extracted from global texture (3 features), Gray-Level Co-occurrence Matrix (GLCM: 9 features), Gray-Level Run-Length Matrix (GLRLM: 13 features), and Neighborhood Gray-Tone Difference Matrix (NGTDM: 5) with 26-voxel connectivity (details of all calculated texture features is presented in the Supplementary Materials, Table S1) using the toolbox developed by Vallière, M. et al. [32]. Prior to TA, normalization and quantization were used to equalize histograms and standardize the intensity range to 64 [32]. quantization were used to equalize histograms and standardize the intensity range to 64 [32]. The chi-square test of independence was used to select the top ten texture features with the highest feature importance. These top ten features were fed into a 5-layered artificial neural network (ANN) classifier model to differentiate between PDAC and non-PDAC masses. We wanted to investigate the role of an ANN, the most basic form of deep neural network, in classification using IVIM-DKI parametric maps because it has shown potential in screening pancreatic cancer [33]. Bayesian optimization was used to optimize the ANN hyperparameters such as layer output size, activation function, lambda, layer weights, and bias initializer. The Bayesian optimizer tuned the hyperparameters to minimize the cross-validation error of the ANN. Following optimization, the output sizes of each hidden layer were set to (221, 16, 3, 2, 2) using the activation function 'tanh', the weights and bias of each layer were started with one using the 'he' initializer, and lambda was set to 0.0061374 for regularization of the hyperparameter tuning. These optimized hyperparameters were fed into an ANN architecture that consisted of an input layer containing predictor variables and class labels, 'tanh' as the activation function for each hidden layer, and a final layer with 'softmax' as the activation function and an output layer that returned the predicted class label. MATLAB was used for the TA and machine learning-based classification (version 9.9, The MathWorks, Inc., Natick, MA, USA).

Statistical Analysis
The patients' mean age and gender distribution were compared using a two-sample t-test and chi-square test. The ADC and IVIM-DKI parameters (D, D*, f, and k) were tested for normality using the Kolmogorov-Smirnov test. The goodness of fit in the standard and IDTV model was assessed using mean adjusted R 2 in tumor ROI, with a high value of adjusted R 2 indicating that the model adequately reflects the data. The coefficient of variation (CV: standard deviation of ROI × 100/mean of ROI) of the IVIM-DKI parameters estimated using standard and IDTV models in pancreatic masses was measured to evaluate the precision of model estimations, where low CV indicates precise estimation of parameters by a model. The Wilcoxon signed-rank test was used to determine significant The chi-square test of independence was used to select the top ten texture features with the highest feature importance. These top ten features were fed into a 5-layered artificial neural network (ANN) classifier model to differentiate between PDAC and non-PDAC masses. We wanted to investigate the role of an ANN, the most basic form of deep neural network, in classification using IVIM-DKI parametric maps because it has shown potential in screening pancreatic cancer [33]. Bayesian optimization was used to optimize the ANN hyperparameters such as layer output size, activation function, lambda, layer weights, and bias initializer. The Bayesian optimizer tuned the hyperparameters to minimize the crossvalidation error of the ANN. Following optimization, the output sizes of each hidden layer were set to (221, 16, 3, 2, 2) using the activation function 'tanh', the weights and bias of each layer were started with one using the 'he' initializer, and lambda was set to 0.0061374 for regularization of the hyperparameter tuning. These optimized hyperparameters were fed into an ANN architecture that consisted of an input layer containing predictor variables and class labels, 'tanh' as the activation function for each hidden layer, and a final layer with 'softmax' as the activation function and an output layer that returned the predicted class label. MATLAB was used for the TA and machine learning-based classification (version 9.9, The MathWorks, Inc., Natick, MA, USA).

Statistical Analysis
The patients' mean age and gender distribution were compared using a two-sample t-test and chi-square test. The ADC and IVIM-DKI parameters (D, D*, f, and k) were tested for normality using the Kolmogorov-Smirnov test. The goodness of fit in the standard and IDTV model was assessed using mean adjusted R 2 in tumor ROI, with a high value of adjusted R 2 indicating that the model adequately reflects the data. The coefficient of variation (CV: standard deviation of ROI × 100/mean of ROI) of the IVIM-DKI parameters estimated using standard and IDTV models in pancreatic masses was measured to evaluate the precision of model estimations, where low CV indicates precise estimation of parameters by a model. The Wilcoxon signed-rank test was used to determine significant differences between mean adjusted R 2 of the two comparative models. A Kruskal-Wallis test with a Tukey-Kramer multi-comparison test was used to compute comparison between IVIM-DKI parameters for subtypes of pancreatic masses (PDAC, pNET, MFCP, and SPEN). The diagnostic performance of a IDTV model was assessed using receiver operating characteristics (ROC) analysis with cut-off values, area under the ROC curve (AUC), sensitivity, and specificity. Important texture features were computed using the chi-square test of independence. The chi-square test of independence was used to determine the relationship between feature variables (texture features derived from ADC or IVIM-DKI parameters) and target variables (class labels: pancreatic masses). If the p-value was small, the feature was considered as an important feature and contributed to one of the significant characteristics in classification; otherwise, it was disregarded if the feature had a large p-value and was not considered as an important feature. Thus, the top ten texture features were chosen to improve performance and simplify the complexity of the classifier model. Performance metrics were estimated such as accuracy, precision, recall, specificity, accuracy, F1 score, and classification error on test data from each stratified fold of the 5-fold cross-validation (training data: n = 38; testing data: n = 10) and this cross-validation was repeated 100 times. A test set from each cross-validation was used to evaluate the performance of the ANN in the classification of PDAC vs. non-PDAC. This 5-fold cross-validation was repeated 100 times; thus, the result was 100 different performance metrics. This validation technique ensures that the data of each patient appeared in the test set a total of 100 times. All statistical analyses were performed using MATLAB (version 9.9, The MathWorks, Inc., Natick, MA, USA).

Model Performance in Pancreatic Masses
Two representative images of IVIM-DKI at high b-value = 2000 s/mm 2 , ADC, and IVIM-DKI parameters maps of a 72-year-old male patient with PDAC are presented in Figure 3a,b,d-k and of a 53-year-old male patient with PDAC in Figure 4a-j and were qualitatively compared. The PDAC-affected region appeared to be hypointense in D* (Figures 3g and 4f) and f maps (Figures 3i and 4h)     The goodness of fit was significantly higher (p < 0.001) for the IDTV model with R 2 = 0.95 ± 0.04 compared to the standard model with R 2 = 0.78 ± 0.16. Figure 3c shows the fitting of the IVIM-DKI signal from one representative patient's data estimated using the standard and IDTV models. Overall, the IDTV model provided precise estimation of IVIM-DKI parameters with lower CV by 22-64% compared to standard model, except for D* in pancreatic masses (the results are presented in the Supplementary Material, Table S2).

Quantitative Comparison between Subtypes of Pancreatic Masses
Figure 5a-e shows the boxplots of the ADC and IVIM-DKI parameters' distribution in PDAC, SPEN, MFCP, and pNET. ADC, D, D*, and k values did not exhibit any significant differences (p > 0.05); however, f values did demonstrate significant differences between pancreatic masses as shown in Table 3. Multiple comparisons between pancreatic masses revealed that the mean f (p < 0.05) was considerably higher in the pNET as compared to the PDAC (the results of multiple comparisons can be seen in Supplementary The goodness of fit was significantly higher (p < 0.001) for the IDTV model with R 2 = 0.95 ± 0.04 compared to the standard model with R 2 = 0.78 ± 0.16. Figure 3c shows the fitting of the IVIM-DKI signal from one representative patient's data estimated using the standard and IDTV models. Overall, the IDTV model provided precise estimation of IVIM-DKI parameters with lower CV by 22-64% compared to standard model, except for D* in pancreatic masses (the results are presented in the Supplementary Material, Table S2).  Table 3. Multiple comparisons between pancreatic masses revealed that the mean f (p < 0.05) was considerably higher in the pNET as compared to the PDAC (the results of multiple comparisons can be seen in Supplementary Material, Table S3)

Differential Diagnosis of Pancreatic Masses Using ROC Analysis
Cut-off values from the ADC and IVIM-DKI parameters were obtained using ROC analysis to differentiate the pancreatic masses such as PDAC, SPEN, MFCP, and pNET, as

Differential Diagnosis of Pancreatic Masses Using ROC Analysis
Cut-off values from the ADC and IVIM-DKI parameters were obtained using ROC analysis to differentiate the pancreatic masses such as PDAC, SPEN, MFCP, and pNET, as shown in Table 4. For PDAC vs. MFCP, only the f parameter showed a high AUC of 0.77 with the cut-off value at 0.20496 having an accuracy, sensitivity, and specificity of 77%, 83%, and 76%, respectively. D* and f showed high AUCs of 0.73 with cut-off 63.73 × 10 −3 mm 2 /s and 0.20424 having an accuracy of 88% and 76%, sensitivity of 63% and 75%, and specificity of 96% and 76%, respectively, to differentiate PDAC from pNET. For pNET vs. MFCP, ADC showed a high AUC of 0.79 with the cut-off 1.58 × 10 −3 mm 2 /s having an accuracy of 79%, sensitivity of 83%, and specificity of 75%. In addition, the D parameter showed a high AUC of 0.76 with the cut-off value at 1.37 × 10 −3 mm 2 /s having an accuracy, sensitivity, and specificity of 79%, 83%, and 75%, respectively.

Multi-Parametric Texture Analysis and Machine Learning-Based Classification of Pancreatic Masses
A total of 30 texture features were calculated for individual parameters, ADC, D, D*, f, and k; thus, 150 textural features (30 × 5) were calculated for each mass. The features of all IVIM-DKI parameters were combined to form a 120 (30 features from each parameter, 30 × 4) texture-feature set. Similarly, ADC and IVIM-DKI parametric maps' texture features were combined to create a total of 150 feature sets. These features were then reduced by selecting only high importance scores which were calculated using the Chi-square test and only the top ten features were selected to reduce the model complexity for PDAC vs. non-PDAC pancreatic masses.
In PDAC vs. non-PDAC, an ANN with all the features of the combined IVIM-DKI parameters showed the highest accuracy of 90.5%, and AUC of 0.92 as shown in Table 5. Classification performance after feature reduction using chi-square decreased by 3-8.5% for combined texture features from IVIM-DKI parameters and individual IVIM-DKI parameters. The top 10 texture features from f showed the highest accuracy of 84.9% and an AUC of 0.85. The combined IVIM-DKI parameters and IVIM-DKI parameters with ADC also performed better with an accuracy of 84.3% and AUC of 0.84. The selected top texture features from the combined IVIM-DKI parameters and IVIM-DKI parameters with ADC were from f, D*, and k. Additionally, decision tree and ensemble were used to characterize PDAC vs. non-PDAC. As shown in Supplementary Table S4, the decision tree performed poorly for PDAC vs. non-PDAC as compared to ensemble and ANN. However, the decision tree showed that the combined texture features of IVIM-DKI parameters after feature reduction showed the highest accuracy and an AUC of 74% and 0.8, respectively. Even the ANN showed a high accuracy and AUC of 84.3% and 0.84 with the combined texture features of IVIM-DKI parameters after feature reduction. Meanwhile, ensemble showed that textural features of the D parameter after feature reduction showed the highest accuracy and an AUC of 91.8% and 0.97, respectively. Commonly selected texture features from multi-parameters were from GLCM (f9), GLRLM (f12, f13, f16, and f21), and NGTDM (f29).

Discussion
MRI can aid in the differential diagnosis of pancreatic masses which is otherwise a very challenging task. With the assumption of a non-Gaussian distribution of water movement, IVIM-DKI can non-invasively measure diffusion and perfusion parameters to capture the tumor heterogeneity. The efficiency of the novel IVIM-DKI model with a total variation penalty function for characterizing pancreatic mass subtypes was evaluated. f was observed to be useful in distinguishing pancreatic masses. In the differential diagnosis of pancreatic masses, f was found to be a crucial indicator in effectively differentiating PDAC from MFCP and both D* and f for differentiating PDAC from pNET. Further, ADC and D could distinguish pNET from MFCP with good accuracy and AUC. In addition, IVIM-DKI parameters were subjected to textural analysis with textural feature reduction and given to the 5-layer artificial neural network for classification of pancreatic masses. PDAC was distinguished from non-PDAC by texture features of f and combined IVIM-DKI parameters with good accuracy and AUC. IVIM-DKI with the IDTV model was able to remove noise from the parametric maps and produce precise parameter values with low CV. These results demonstrate that the IVIM-DKI model with a spatially constrained optimization method may provide true underlying texture features for successful characterization of pancreatic masses.
A new approach was used to precisely measure IVIM-DKI parameter values to differentiate between pancreatic masses. This IDTV model has been effectively applied in characterization of prostate lesions [16] and patients with lymphoma [34]. The TV penalty function is an adaptive regularization approach that iteratively balances the regularization weight while preserving lesion edge [17,18]. The IDTV model has demonstrated good description of the underlying data in tumor tissues and reduced fitting residuals when estimating the parameters, which was consistent with the findings of Baidya Kayal et al. (2017) [18], who also demonstrated improved parameter estimations using the IVIM model with TV in osteosarcoma [18]. Repeatability was not assessed in this study because we did not obtain repeat scans; however, previous research has shown significant inter-scan agreement and no variability in parameter estimation in healthy tissue volume across time points, proving the repeatability and reproducibility of this new approach in clinical scenarios [35]. Additionally, parametric maps assessed using the IDTV model were found to be less noisy than the standard model [13,16,18]. However, the real clinical impact of these parametric maps estimated using the IDTV model would need multi-centric and multi-rater assessment.
Although the efficacy of IVIM-DKI in pancreatic mass characterization has yet to be completely examined, it provides a quantitative evaluation of perfusion and heterogeneity in the tumor. Mean f was able to characterize PDAC from pNET, where f was high in the latter, which might be attributed to the large amount of desmoplastic stroma and poor tumor vascularity in PDAC [36][37][38][39]. Meanwhile, pNET has a well-defined capillary network and is a hypervascular tumor compared to PDAC [36,37]. ADC, D, D*, and k were unable to distinguish between pancreatic masses, which was similarly demonstrated in previous studies [30,36,37,40,41]. This might be due to the differentiating feature of PDAC being tissue perfusion rather than diffusion [30]. Even when using the IDTV model, D* exhibited a high estimation error, indicating that it cannot be employed as a stand-alone biomarker due to its high sensitivity to signal noise and limited repeatability [36,42].
Additionally, this study estimated cut-off values for further subtyping of pancreatic masses. ADC and D showed good accuracy and AUC for characterizing pNET vs. MFCP. The presence of granulation and fibrosis in MFCP [43] may have contributed to the good diagnostic performance of ADC and D for differentiating pNET from MFCP. The D* and f parameters demonstrated a good degree of accuracy and AUC in discriminating PDAC from pNET; f was also useful for discriminating PDAC from MFCP. However, Klauss et al. (2011) [41] showed that the f parameter had a higher AUC in differentiating between PDAC and MFCP than this study (using a IVIM-DKI model with regularization). The former study used a fixed value of D* derived from the healthy subjects and used it in the model to estimate parameter f; this approach was adopted by Klauss et al. to avoid unreliable parameter estimation (f value) by the simultaneous fitting of all four parameters. However, for patient data, this technique may result in a variation in the estimation of f, as D* has been fixed based on a healthy population. In comparison to their model, the IDTV model does not fix or have prior assumptions for any parameter and estimates all the four parameters concurrently which was also done in previous studies using the same model in prostate cancer [16] and lymphoma [34].
Chronic pancreatitis of the mass-forming type has an imaging appearance like PDAC, leading to misdiagnosis on cross-sectional imaging. On CT perfusion imaging (CTP), blood volume (BV) and blood flow (BF) have been shown to differentiate PDAC from MFCP [44,45]. Perfusion fraction (f) has a significant correlation with CTP parameters such as BV and BF [46,47] and microvessel parameters which are an indicator of tumor angiogenesis in PDAC [47]. This further supports that the IVIM-DKI parameters could be an alternative in patients with compromised kidney function where administration of iodinated or gadolinium-based contrast agents is contraindicated. In comparison to earlier studies, only parameter f demonstrated a greater accuracy AUC for discriminating PDAC from MFCP [36]. While it is widely known that differentiating PDAC from MFCP can be problematic because of the overlapping imaging characteristics, perfusion fraction can be effective in accomplishing this task.
In assessing pancreatic masses, whole-volume TA was also carried out for PDAC vs. other masses. TA can offer information on non-invasively detecting small variations in the intensity of a lesion location and can yield additional information on microscopic tissue heterogeneity [40]. This study used texture features of individual and combined IVIM-DKI parameters and the combination of IVIM-DKI parameters and ADC texture features. PDAC was successfully distinguished from other pancreatic tumors utilizing textural features from f and integrated IVIM-DKI parameters. ADC contributed minimal value to this classification performance for detecting PDAC. The TA of individual IVIM-DKI parameters such as f outperformed other parameters after feature reduction in the classification. Following feature reduction, texture features retrieved from integrated IVIM-DKI parameters largely comprised D*, f, and k features and features of GLCM, GLRLM, and NGTDM methods. The texture features from k performed fairly with good diagnostic performance in detecting PDAC; however, the advantages of the kurtosis parameter have not been widely investigated in the literature [38]. Only a few studies have looked at IVIM parameter texture analysis [23], and none have used higher-order texture analysis of IVIM-DKI parameters such as the NGTDM method. An ANN has been used in differentiation of pancreatic masses using dynamic contrast-enhanced MRI but not with IVIM-DKI parameters [48].

Limitation and Future Scope
This study has a few limitations including, first, the study cohort was small, and the results from analysis might be biased due to it being a single center study and thus conclusion can be considered as preliminary evidence requiring wider validation. Further, this study needs to be performed with large multicentric data analysis. Second, in this cohort, age and gender were not matched; however, no correlation between age or gender distribution with ADC or IVIM-DKI parametric values were observed. Third, the acquisition time to acquire IVIM-DKI was 16 min, which might cause motion artifacts; hence, future investigations will be required to optimize the b-values for a shorter acquisition duration. Fourth, the correlation between IVIM-DKI parameters and histopathological markers were not assessed in the current study. Fifth, the IDTV model was not compared with Bayesian techniques since it was outside the scope of this study. Lastly, the texture analysis comparing all possible subtypes of pancreatic masses was not performed because the sample size was insufficient to provide conclusive evidence for subtypes analysis; however, data augmentation techniques were used for this limitation, but results were found to be over-biased and inconclusive due to very small dataset for MFCP and SPEN.
A larger dataset of rare pancreatic cancers such as SPEN and MFCP needs to be acquired so that this can help us in understanding its contribution in characterization of pancreatic masses. Machine learning-based texture analysis may benefit not only in the accurate diagnosis and management of cancer, but also in the virtual reality-based surgical procedure planning [49,50]. This method has the potential to aid surgeons in the navigation of hidden structures with few incisions, such as the precise delineation of dissection planes, such as the pancreatic tumor parenchyma. Automatic segmentation of pancreatic masses can be achieved using deep learning with multi-centric large datasets for accurate delineation of tumors.

Conclusions
In conclusion, this study demonstrated that the IVIM-DKI model with the spatial regularization optimization approach has a strong potential to be used as an additional paradigm for differential diagnosis of pancreatic masses. Perfusion measures such as D* and f were found to out-perform with accuracies of 88% and 77%, respectively in discriminating pancreatic adenocarcinoma from pancreatic neuroendocrine tumor and mass forming-chronic pancreatitis than the diffusion measures (ADC and D). Non-invasive characterization of pancreatic masses can also largely benefit from whole-volumetric texture analysis of f and combined texture features from IVIM-DKI parameters paired with machine learning-based classification with accuracies of 84.9% and 84.3%, respectively.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/bioengineering10010083/s1, Figure S1: 1 to 14: DWI trace images of a representative patient from IVIM-DKI images at 14 b-values from 0 to 2500 s/mm 2 (0, 25, 50, 75, 100, 150, 200, 500, 800, 1000, 1250, 1500, 2000, 2500 s/mm 2 ); Figure S2: Representative image of a 28-year-old female patient with solid papillary epithelial neoplasm (SPEN) with qualitative comparison between parameter maps (D, D*, f, and k maps) estimated using standard and IDTV models; Figure S3: Representative image of a 26-year-old female patient with solid papillary epithelial neoplasm (SPEN) with qualitative comparison between parameter maps (D, D*, f, and k maps) estimated using standard and IDTV models; Figure S4: Representative image of a 63-year-old male patient with mass-forming chronic pancreatitis (MFCP) with qualitative comparison between parameter maps (D, D*, f, and k maps) estimated using standard and IDTV models; Figure S5: Representative image of a 60-year-old male patient with mass-forming chronic pancreatitis (MFCP) with qualitative comparison between parameter maps (D, D*, f, and k maps) estimated using standard and IDTV models; Figure S6: Representative image of a 43-year-old male patient with pancreatic neuroendocrine tumor (PNET) with qualitative comparison between parameter maps (D, D*, f, and k maps) estimated using standard and IDTV models; Figure S7: Representative image of a 41-year-old male patient with pancreatic neuroendocrine tumor (PNET) with qualitative comparison between parameter maps (D, D*, f, and k maps) estimated using standard and IDTV models; Figure S8: Representative image of a 33-year-old male patient with pancreatic ductal adenocarcinomas (PDAC) with qualitative comparison between parameter maps (D, D*, f, and k maps) estimated using standard and IDTV models; Table S1: Texture features extracted from Global, GLCM, GLRLM, and NGDZM; Table S2: Coefficient of variation (CV: standard deviation of ROI × 100/mean of ROI) of IVIM-DKI parameters estimated using standard and IDTV models in pancreatic masses; Table S3: p-values obtained from Tukey-Kramer multiple comparison test between ADC and IVIM-DKI parameters values estimated using IDTV model in pancreatic masses; Table S4: Machine learning-based classification performance of all texture features (30 features) and top ten features selected from ADC, combined, individual IVIM-DKI parameters, and ADC combined with IVIM-DKI parameters to classify PDAC from non-PDAC pancreatic masses using decision tree and ensemble.  Informed Consent Statement: Informed consent was obtained from all the subjects involved in the study. Written informed consent was obtained from the patient(s) to publish this paper.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author (A.M.) upon reasonable request.