Gamma distribution model of diffusion MRI for the differentiation of primary central nerve system lymphomas and glioblastomas

The preoperative imaging-based differentiation of primary central nervous system lymphomas (PCNSLs) and glioblastomas (GBs) is of high importance since the therapeutic strategies differ substantially between these tumors. In this study, we investigate whether the gamma distribution (GD) model is useful in this differentiation of PNCSLs and GBs. Twenty-seven patients with PCNSLs and 57 patients with GBs were imaged with diffusion-weighted imaging using 13 b-values ranging from 0 to 1000 sec/mm2. The shape parameter (κ) and scale parameter (θ) were obtained with the GD model. Fractions of three different areas under the probability density function curve (f1, f2, f3) were defined as follows: f1, diffusion coefficient (D) <1.0×10−3 mm2/sec; f2, D >1.0×10−3 and <3.0×10−3 mm2/sec; f3, D >3.0 × 10−3 mm2/sec. The GD model-derived parameters were compared between PCNSLs and GBs. Receiver operating characteristic (ROC) curve analyses were performed to assess diagnostic performance. The correlations with intravoxel incoherent motion (IVIM)-derived parameters were evaluated. The PCNSL group's κ (2.26 ± 1.00) was significantly smaller than the GB group's (3.62 ± 2.01, p = 0.0004). The PCNSL group's f1 (0.542 ± 0.107) was significantly larger than the GB group's (0.348 ± 0.132, p<0.0001). The PCNSL group's f2 (0.372 ± 0.098) was significantly smaller than the GB group's (0.508 ± 0.127, p<0.0001). The PCNSL group's f3 (0.086 ± 0.043) was significantly smaller than the GB group's (0.144 ± 0.062, p<0.0001). The combination of κ, f1, and f3 showed excellent diagnostic performance (area under the curve, 0.909). The f1 had an almost perfect inverse correlation with D. The f2 and f3 had very strong positive correlations with D and f, respectively. The GD model is useful for the differentiation of GBs and PCNSLs.


Introduction
The preoperative imaging-based differentiation of primary central nervous system lymphomas (PCNSLs) and glioblastomas (GBs) is of high importance since the therapeutic strategies differ substantially between these tumors [1,2]. The treatment of GBs is based on the maximal possible safe surgical resection together with postoperative chemoradiation therapy [1] whereas PCNSLs require a biopsy for histological confirmation followed by chemoradiation therapy [2]. In typical cases, the differentiation of these tumors by conventional MRI is not always difficult since PCNSLs shows homogenous contrast enhancing lesions while GBs show irregular and heterogenous ring enhancing lesion with necrosis. However, it is frequently difficult to differentiate these tumors especially when they demonstrate atypical imaging features.
Many mathematical models have been proposed for the analysis of diffusion MRI. The mono-exponential model describes the Brownian motion of water molecules by calculating the apparent diffusion coefficient (ADC) based on the Gaussian distribution of diffusion displacement [3]. The bi-exponential intravoxel incoherent motion (IVIM) model aims to separate the true water diffusion and the capillary perfusion by using multiple low b-values [10,11]. Diffusion kurtosis imaging (DKI) is an approach used to characterize non-Gaussian water diffusion, which estimates kurtosis metrics [12].
It has been reported that all of these approaches are useful in differentiating GBs and PCNSLs [3,13,14], but all have possible limitations. The mono-exponential model may not precisely reflect the reality of diffusion behavior in heterogenous biological tissues, since this model assumes a Gaussian distribution. The bi-exponential model could be influenced by an uncertainty of the estimated perfusion, since signal measurements at low b-values are susceptible to measurement errors [15][16][17][18]. The DKI model is limited by the unclear biological interpretation of the kurtosis parameters [18][19][20][21].
As one of the non-Gaussian distribution models, a statistical model based on the gamma distribution (GD) has been proposed for diffusion MRI analyses [22]. The GD model is a twoparameter family of continuous probability distribution parametrized in terms of the shape parameters kappa (κ) and the scale parameter theta (θ), and this model assumes that the diffusion coefficient (D) is distributed continuously within a voxel. The GD model allows us to estimate fractions of a tissue type based on the concept that the area fractions for D <1.0 × 10 −3 mm 2 /sec, D = 1.0 × 10 −3 to 3.0 ×10 −3 mm 2 /sec, and D >3.0 ×10 −3 mm 2 /sec are attributed to intracellular, extracellular extravascular, and intravascular spaces, respectively [18,22,23]. Based on these fractions, we may be able to estimate histopathological conditions of neoplasms or organs.
The GD model has been used to assess prostate cancers [22][23][24], breast cancers [18], and renal function [25]. The GD model was also used to assess cerebral ischemic stroke in rat brains, and it was showed that this model exhibited a better performance than the conventional mono-exponential model and allowed for a significantly enhanced visualization of ischemic lesions [26]. To the best of our knowledge, its application to brain tumors has never been reported. We conducted the present study to determine whether the GD model is useful in the differentiation of PCNSLs and GBs.

Materials and methods
This retrospective study was approved by the Institutional Review Board of Kyushu University Hospital (no. 2019-447), and the requirement for informed consent was waived.

Patients
The DWI protocol with multiple b-values has been a part of our routine preoperative MRI examination for patients with brain tumors since January 2013. The patient inclusion criteria for this study were: (1) The DWI with multiple b-values was conducted preoperatively for the patient during the period from January 2013 to August 2019; and (2) The patient subsequently underwent a surgical resection or biopsy within 1 month of the DWI with multiple b-values, and the histopathological diagnosis of PCNSL or GB was made. A total of 89 patients met these criteria. The exclusion criteria were as follows: (1) no distinct contrast enhancement observed in the lesion (n = 3); and (2) difficulty in the analysis of images due to severe artifacts (n = 2). Thus, a total of 84 patients including 27

Image analysis
The mono-exponential model was computed using all of the above-listed b-values according to the following equation: where Sb is the signal intensity for each b-value and S0 is the signal intensity at a b-value of zero.
In the bi-exponential model, the signal decay was estimated by the following the equation: where D � is the pseudo-diffusion coefficient, and the f is the volume fraction within a voxel of water flowing in perfused capillaries. The GD model is represented by ρ(D) and is given by: where κ describes the shape parameter and θ describes the scale parameter. When the distribution of D follows this equation, the signal intensity on DWI is given by: Three different areas under the probability density function (PDF) curve were defined as follows: f1, the fraction of D <1.0×10 −3 mm 2 /sec; f2, the fraction of 1.0×10 −3 to 3.0×10 −3 mm 2 / sec; f3, the fraction of D >3.0×10 −3 mm 2 /sec. The f1 value is attributed to the intracellular component; the f2 is attributed to the extracellular extravascular component, and the f3 is attributed to the intravascular component [18,22,23].
The DWI data in the digital imaging and communications in medicine (DICOM) format were transferred to a personal computer and fit to the GD model, and then the κ and θ values were estimated using the Image J software program (ver. 1.52p; U.S. National Institutes of Health, Bethesda, MD) and self-built plug-ins. After the export of the x-and y-coordinates and the κ and θ of each pixel within the region of interest (ROI), the f1, f2, and f3 values of each pixel were calculated using Microsoft Excel ver. 16.16.14.

ROI placement
The matrix sizes of the postcontrast T1-weighted images were adjusted to the same size as those of the DWI using the ImageJ function to match the geometric information of these images. ROIs were placed to delineate the enhancing lesion on the single slice that had the maximum area. On the size-adjusted postcontrast T1-weighted images, areas showing contrast enhancement were manually segmented by a neuroradiologist with 19 years of experience (O. T.) (Fig 1). The areas with necrosis, cystic lesion, hemorrhage, or obvious artifacts were carefully excluded from the ROI.
The ROIs were copied from the postcontrast T1-weighted images and pasted to the DWI. Fine manual adjustments were made when there were locational mismatches due to image distortion or the patient's motion, etc. The ROIs were also placed on the peritumoral non-contrast-enhancing T2-hyperintense areas to evaluate whether there were differences in histological features including tumor infiltration or increased vascularity in the peritumoral areas between PCNSLs and GBs. In addition, the ROIs were placed on the contralateral normal-appearing white matter. The ROIs for the peritumoral non-contrast-enhancing T2-hyperintense areas and contralateral normal-appearing white matter were measured on the image obtained with the b-value of 0 sec/mm 2 image. The same ROIs were used for all DWI analyses.

Statistical analyses
The GD model-derived and IVIM-derived parameters were compared between the PCNSLs and GBs with the Mann-Whitney U-test. A receiver operating characteristic (ROC) curve analysis was performed to assess the diagnostic performance of the parameters in the differentiation of PCNSLs and GBs. The area under the curve (AUC) was calculated, and then the sensitivity and specificity were obtained. The optimal cutoff point was determined by Youden's method [28]. The diagnostic performance was considered excellent for AUC values between 0.9 and 1.0, good for AUC values between 0.8 and 0.9, fair for AUC values between 0.7 and 0.8, poor for AUC values between 0.6 and 0.7, and failed for AUC values between 0.5 and 0.6 [29].
To determine whether the combination of multiple parameters for both the GD model and the IVIM model increases the diagnostic performance, we first performed a stepwise analysis to select the explanatory variables for a multiple regression model from a group of candidate variables by going through a series of automated steps. A forward-selection rule was applied in which the analysis started with no explanatory variables and then added variables, one by one, based on which variable was the most statistically significant, until there were no remaining statistically significant variables [30,31]. We then performed a binomial logistic regression analysis to examine the AUCs of the combinations of the selected parameters. Two independent AUCs were compared using the method of Delong et al. [32]. The correlations among the parameters were assessed with Pearson' correlation. Statistical analyses were performed with Prism 5.0 (GraphPad Software, San Diego, CA), MedCalc 19.1 (Broekstraat, Mariakerke, Belgium), and JMP Pro 14.0 (SAS Institute, Cary, NC). P-values <0.05 were considered significant.

Comparisons of the parameters between the PCNSL and GB groups
The detailed information for the parameters in the gadolinium enhancing lesion, peritumoral T2-hyperintense areas without contrast enhancement, and normal appearing white matter is summarized in Table 1.
The results of our comparisons of the GD model-derived parameters between the PCNSLs and GBs in the gadolinium-enhancing lesions are shown in Fig 2. In the gadolinium-enhancing lesions, the κ was significantly smaller in the PCNSL group (2.26 ± 1.00) than in the GB group (3.62 ± 2.01, p = 0.0004), the f1 was significantly larger in the PCNSL group (0.542 ± 0.107) than in the GB group (0.348 ± 0.132, p<0.0001), the f2 was significantly smaller in the PCNSL group (0.372 ± 0.098) than in the GB group (0.508 ± 0.127, p<0.0001), and the f3 was significantly smaller in the PCNSL group (0.086 ± 0.043) than in the GB group (0.144 ± 0.062, p<0.0001), while the θ was not significantly different between the groups.
In the peritumoral T2-hyperintense areas without contrast enhancement, no significant differences were found between the PCNSL and GB groups for any of the GD model derived parameters.

Diagnostic performance of the single and combined parameters
The ROC graphs and diagnostic performance parameters are shown in Fig 5 and Table 2. In the single-parameter analysis regarding the differential diagnosis of GBs and PCNSLs, the  In the combined-parameters analysis, the stepwise procedure selected κ, f1, and f3 for the GD model, and the D and f for the IVIM model. The combination of κ, f1, and f3 revealed excellent diagnostic performance with the AUC of 0.909, sensitivity of 84.2%, and specificity of 88.9%. This combination increased the diagnostic performance of κ (p = 0.0016), and f3 (p = 0.0075), although it did not improve the performance of f1 (p = 0.1950). The AUC of this combination (0.909) was higher than that of ADC (0.879); however, there was no significant difference between them (p = 0.2152).
The combination of D and f showed good diagnostic performance with the AUC of 0.884, 82.5% sensitivity, and 81.5% specificity. This combination improved the diagnostic performance of f (p = 0.0077), although it did not improve the performance of D (p = 0.5276).
Among all of the single and combined parameters, the combination of κ, f1, and f3 showed the highest AUC; however, no significant differences were detected between this combination and the ADC (p = 0.2152) or the combination of D and f (p = 0.2207).

Discussion
The results of our analyses revealed that in gadolinium-enhancing lesions, the κ was significantly smaller in the PCNSL group than in the GB group. The θ was not different between the groups. The f1 was larger, the f2 was smaller, and the f3 was lower in the PCNSLs than in the GBs. The low κ values observed in the PCNSLs indicated that the PDF curve had a rightskewed distribution, which meant that the PDF has its peak in the lower D area, and thus the fraction of lower D was larger. Since the θ values were not significantly different between the PCNSL and GB groups, it was likely that the lower κ values might result in the lower ADC and D values and the higher f1 values observed in the PCNSLs compared to the GBs. These findings are in accordance with studies that examined the mono-exponential model, in which PCNSLs showed lower ADC values relative to GBs [3][4][5].
The θ is a scale parameter and may thus reflect the heterogeneity of a biological tissue. We expected that the θ values would be larger in GBs than in PCNSLs since GBs are histologically characterized by intratumoral tissue heterogeneity whereas PCNSLs are characterized by the dense and homogenous distribution of tumor cells; however, no significant difference in the θ values was observed between the groups. The θ values showed large standard deviations in both the PCNSLs and the GBs, indicating that this value could vary widely even in the same type of tumor. The same trend was observed in a study of breast tumors in which the θ values were not significantly different between the different types [18]. The utility of this parameter should be further evaluated in larger populations.
It seems that the higher f1 and lower f2 in the PCNSLs and the lower f1 and higher f2 in the GBs well reflected the histological features of the respective tumors. Histologically, PCNSLs are characterized by high cell density at the expense of reduced available extracellular space, and necrosis is not a common feature of this tumor. GBs can show locally high cell density, but the overall cell density can be lowered depending on the fraction of microscopic necrosis or hemorrhage. Our present findings are consistent with a study that reported that the ADC was lower and the cell density was higher in PNCSLs than in high-grade gliomas [3].
The GB group showed larger f3 and f compared to the PCNSL group. This may be attributed to the difference in vascularity of these tumors. Pathologically, neovascularization is a key feature of GB while it is not prominent in PCNSL [33,34]. Our results are consistent with those from previous studies using dynamic susceptibility contrast perfusion-weighted imaging and arterial spin labeling imaging [9,35].
With respect to the diagnostic performance, the ADC, f1, and D showed comparable AUCs in the present study, indicating that all three of these parameters are useful in the differentiation of PNCSLs and GBs. The reason for the slightly higher AUC observed with the ADC could be the effect of perfusion on ADC measurements. In hyperperfused tissues, ADC will be affected by the perfusion effect and overestimated compared to D; however, since both f1 and D are parameters without a perfusion effect in theory, an overestimation caused by perfusion should not be observed in these values. Therefore, in hypervascular tumors such as GBs, the ADC should be larger than D. On the other hand, in hypovascular tumors such as PCNSLs, this difference between ADC and D should smaller. This means that the difference between ADC and D would be larger in GBs than in PCNSLs. Therefore, ADC could show higher diagnostic performance in the discrimination of these two tumors than D. In fact, the difference between the ADC and D values was greater in the GBs (0.100 × 10 −3 mm 2 /sec) than in the PCNSLs (0.078 × 10 −3 mm 2 /sec), which was most likely due to the higher perfusion effect on the ADC in GBs than in PCNSLs. Nevertheless, the combination of κ, f1, and f3 demonstrated the highest diagnostic performance among all of the single and combined parameters, with the AUC of 0.909. The AUC of this combination tended to be higher than that of ADC although there was no statistically significant difference. Whether the combination of parameters of the GD model has an additive value should be evaluated in a larger population, since we did not observe statistical significance in all of our comparisons.
We found the correlations between the GD model-derived and IVIM-derived parameters, particularly between the f1 and D, the f2 and D, and the f3 and f. The almost perfect negative correlation observed between the f1 and D may indicate that these two parameters contain virtually identical information. The positive correlation between f2 and D suggests that the increased extracellular space like that taken up by microscopic necrosis might result in the higher f2. The positive correlation between f3 and f indicates that both of these parameters well reflected tissue perfusion despite the different analysis methods used. The negative correlation between f1 and f2 was likely due to the complementary relationship between these two parameters. In general, intravascular space (≒ f3) is smaller compared to intracellular (≒ f1) and extracellular extravascular space (≒ f2). In fact, the f3-values were much smaller than the f1-and f2-values in both PCNSLs and GBs in the present study. Therefore, the increase in f1 would result in the decrease in f2, and vice versa. Although the GD-derived and IVIM-derived parameters provide similar information, the strength of the GD model-derived parameters is that all fraction values (f1, f2, f3) are expressed as fractions or percentages, which allows us to well characterize tumors from histological viewpoint. The IVIM-derived f-value is also expressed in a percentage or fraction; however, the IVIM analysis is not able to provide the fraction values for intracellular and extracellular-extravascular spaces. In this sense, the IVIM method is not a perfect method for the histological characterization of tumors.
In the T2-hyperintense lesions without contrast enhancement, no significant differences were observed between the PCNSL and GB groups for any parameters. There have been several studies that showed increased rCBV on DSC-perfusion imaging in peritumoral noncontrast-enhancing T2-hyperintense areas of GBs [36,37]. The results of these studies indicated that the peritumoral areas of GB include not only vasogenic edema but also tumor cells infiltrating surrounding brain parenchyma; however, our study did not reveal any significant differences in the GD model-based parameters for peritumoral noncontrast-enhancing T2-hyperintense areas between PCNSLs and GBs. The f2 values in the noncontrast-enhancing T2-hyperintense areas were higher in both types of tumor compared to those in the contrastenhancing areas and normal appearing white matter. We assume that the high f2 values in the noncontrast-enhancing T2-hyperintense areas are likely to reflect mostly perifocal vasogenic edema rather than tumor infiltration outside the enhancing lesion. Our result is consistent with the previous DWI study in which ADC could not be used to differentiate edema with infiltration of tumor cells from vasogenic edema in high-grade gliomas and PCNSLs [38].
In the normal-appearing white matter, the GB group showed larger f1, smaller f2, and larger f3 than the PCNSL group although these differences were small. This was unexpected, and the reasons for the differences remain unclear; however, since GBs frequently show extensive infiltration into the surrounding brain tissue, which is a fundamental feature of diffuse glioma, it is no wonder that the increased cell density and perfusion were observed in the normal-appearing white matter.
This study has several limitations. The number of patients was relatively small (n = 84)especially the number of patients with PCNSL (n = 27). The only one person performed the ROI placements on a single slice, and not whole tumor volume was evaluated. The ROI placements on the gadolinium-enhancing lesions were occasionally difficult, particularly when the lesions showed irregular and thin ring-like enhancement. Although the best effort was made to include only enhancing lesions, it is possible that necrosis in tumors was included, and this could have affected the analyses. In addition, the selection of b-values has not yet been optimized. Prior studies of the GD model used the maximum b-values ranging from 1000 to 3000 sec/mm 2 [18,[22][23][24]. In a study of prostate cancers, Oshio et al. used the similar DWI parameters to ours and the highest b-value of 1000 s/mm 2 , and reported that the good fitting accuracy was observed in both cancerous tissues (R 2 = 0.99226) and normal tissues (R 2 = 0.99842) [22]. Their result indicated that DWI with the highest b-value of 1000 s/mm 2 can be used for GD model analyses; however, since it was reported that the non-monoexponential diffusionrelated signal decay generally becomes more obvious over more extended b-value ranges, the maximum b value of 1000 sec/mm 2 used in the present study might be lower than the optimal value. The optimal b-values and numbers should be elucidated in future studies.

Conclusions
The GD model well described the histological features of PCNSLs and GBs, and its use enabled the significant differentiation of these tumors. The κ, f2, and f3 values were significantly smaller and the f1 values were significantly larger in the PCNSLs than in the GBs. The combination of κ, f1, and f3 showed the highest AUC. The GD model-derived parameters correlated well with the IVIM-derived parameters. The GD model may therefore contribute to the characterization of various brain tumors from the histological viewpoint.