Morphological, Functional and Texture Analysis Magnetic Resonance Imaging Features in the Assessment of Radiotherapy-Induced Xerostomia in Oropharyngeal Cancer

: The aim of this single-center, observational, retrospective study was to investigate magnetic resonance imaging (MRI) biomarkers for the assessment of radiotherapy (RT)-induced xerostomia. Twenty-seven patients who underwent radiation therapy for oropharyngeal cancer were divided into three groups according to the severity of their xerostomia—mild, moderate, and severe—clinically conﬁrmed with the Common Terminology Criteria for Adverse Events (CTCAE). No severe xerostomia was found. Conventional and functional MRI (perfusion-and diffusion-weighted imaging) performed both pre-and post-RT were studied for signal intensity, mean apparent diffusion coefﬁcient (ADC) values, k-trans, and area under the perfusion curves. Contrast-enhanced T1 images and ADC maps were imported into 3D slicer software, and salivary gland volumes were segmented. A total of 107 texture features were derived. T-Student and Wilcoxon signed-rank tests were performed on functional MRI parameters and texture analysis features to identify the differences between pre-and post-RT populations. A p -value < 0.01 was deﬁned as acceptable. Receiver operating characteristic (ROC) curves were plotted for signiﬁcant parameters to discriminate the severity of xerostomia in the pre-RT population. Conventional and functional MRI did not yield statistically signiﬁcant results; on the contrary, ﬁve texture features showed signiﬁcant variation between pre-and post-RT on the ADC maps, of which only informational measure of correlation 1 (IMC 1) was able to discriminate the severity of RT-induced xerostomia in the pre-RT population (area under the curve (AUC) > 0.7). Values lower than the cut-off of − 1.473 × 10 − 11 were associated with moderate xerostomia, enabling the differentiation of mild xerostomia from moderate xerostomia with a 73% sensitivity, 75% speciﬁcity, and 75% diagnostic accuracy. Therefore, the texture feature IMC 1 on the ADC maps allowed the distinction between different degrees of severity of RT-induced xerostomia in the pre-RT population. Accordingly, texture analysis on ADC maps should be considered a useful tool to evaluate salivary gland radiosensitivity and help identify patients at risk of developing more serious xerostomia before radiation therapy is administered.


Introduction
Head and neck cancers represent the sixth most common form of cancer worldwide. More than 90% of head and neck cancers are squamous cell carcinomas of the oral cavity, oropharynx, and larynx mucosal tissues [1]. Radiotherapy (RT) is generally included in the primary oncologic treatment, as it improves clinical and functional outcomes for cancer patients, especially in the case of oropharyngeal cancers (OPC) [2]. In addition to its therapeutic effects, RT can cause injuries to the tissues inside and surrounding the irradiated area, and approximately 90% of patients undergoing radiotherapy for head and neck cancers suffer from clinically relevant xerostomia [3,4].
Xerostomia is defined as dryness of the oral cavity resulting from insufficient saliva secretion [5,6]. RT-induced xerostomia is caused by salivary gland dysfunction resulting from X-ray related tissue damage [5,7,8]. In regard to the diagnosis of xerostomia, the basic tests include the determination of stimulated and unstimulated salivary flow rate, palatal secretion, and parotid secretion [9][10][11]. These measurements constitute the simplest methods of assessing the secretory function of salivary glands. Very low unstimulated and stimulated salivary flow rates are defined as <0.1 mL/min and <0.7 mL/min, respectively [10][11][12]. Such values are confirmatory of xerostomia, whether or not they co-exist with specific symptoms for this condition, such as oral soreness, dry lips, halitosis, decreased or altered sense of taste, recurrent mouth infections, tooth decay and gum disease, and difficulty speaking, eating, or swallowing [13][14][15].
Because of the high-contrast resolution and the ability to study complex anatomical regions without the use of radiation, magnetic resonance imaging (MRI) is considered the most relevant imaging technique for the identification of head and neck lesions [16][17][18][19]. Diffusion-weighted imaging (DWI) is an established diagnostic tool that evaluates the tissue microanatomy by studying the spontaneous molecular diffusion of protons corresponding to the stochastic Brownian motion of water molecules [16,20,21]. The apparent diffusion coefficient (ADC) is calculated from DWI and consists of the quantitative assessment of the impedance of water molecule diffusion within tissues [20]. The ADC is typically reduced in hypercellular tissues and increased in situations where water molecules are free to move [16,17]. The ADC values are calculated automatically and integrated into a parametric map, upon which regions of interest (ROIs) can be traced at a workstation to determine the ADC values of specific portions of a tissue [16][17][18]. Dynamic contrastenhanced, perfusion-weighted imaging (DCE-PWI) is a functional MRI technique that allows one to infer blood perfusion to specific tissues by measuring the changes in tissues over time after the intravenous administration of a contrast agent [22]. DCE-PWI is crucial in the detection of focal lesions since it mainly assesses the vascular permeability [23]. Changes in hemodynamic parameters can precede abnormalities on conventional MRI and can thus be used to help with the diagnosis [22][23][24][25][26].
Radiomics has been rapidly developing over the last few years; it is a hybrid analytical process aimed at determining the correlation between the characteristics of tissues and their corresponding digital images [17,27]. Texture analysis is a form of radiomics, in which macroscopic heterogeneities of tissues can be non-invasively studied to infer information about their microscopic architecture beyond the possibilities of the human eye as if it were a "virtual biopsy" [17,28,29]. Texture analysis is based on the extraction of parameters representing the distribution frequency, intensity, or direction of gray levels within the ROI in order to evaluate the single pixel, its interactions with adjacent pixels, and the distribution of pixels and voxels in the image [17,29].
Some previous studies evaluated RT-induced xerostomia with different MRI techniques, such as DWI [21,[30][31][32] and DCE-PWI [33]. However, neither used both techniques on the same cohort of patients in pre-and post-RT MRI examinations. Furthermore, the texture analysis was mainly performed with CT imaging [34][35][36], and only one paper can be found on ultrasound [37] and conventional T1 MRI sequences [38].
The current study represents the first attempt to evaluate RT-induced xerostomia by using multiparametric MRI techniques, including DWI, DCE-PWI, and texture analysis, carried out in both pre-and post-RT imaging. The aim of this retrospective study was to correlate such MRI techniques with the severity of RT-induced xerostomia clinically confirmed with the Common Terminology Criteria for Adverse Events (CTCAE). The initial population included 180 patients who underwent RT to treat oropharyngeal cancer ( Figure 1). Among them, 128 did not have available MRI data. Of the remaining 52 patients who were studied with MRI, 21 patients had no DWI and DCE-PWI performed or MRI examinations both before and after RT were not carried out. In addition, 4 patients were not clinically assessed or did not develop clinically confirmed xerostomia. Therefore, the final number of patients included in the study was 27. They developed xerostomia after RT and were studied with MRI both before and after RT with DWI and DCE-PWI sequences.

Patients' Differentiations into Groups Based on Clinical Evaluation
The severity of dry mouth was clinically assessed 1 month after the end of RT by a radiotherapist with 10-years' experience using the U.S. National Cancer Institute's CTCAE v4. This score scale has three increasing levels of severity [39]. Accordingly, we decided to divide the patients into three groups:

Image Acquisition and Analysis
MRI examinations were performed via 1.5 T Magnetom aera (Siemens Healthcare, Erlangen, Germany) with a devoted head and neck coil. The MRI acquisition protocol included pre-and post-contrast scans. An axial fat-saturated, echo-planar, imaging-based DWI with two different b-values (b 50-800 s/mm 2 ) was acquired. The apparent diffusion coefficient (ADC) values of the parotid and sub-mandibular glands before and after RT were calculated by positioning three regions of interest (ROI) within the salivary gland parenchyma in three contiguous axial sections ( Figure 2). Time/intensity curve (I/t), area under the curve (AUC), and K(trans) values of the parotid and submandibular glands before and after RT were generated by using IntelliSpace software version 9.0 (Philips, Amsterdam, The Netherlands) from the native DCE-PWI images by drawing a ROI including the largest gland section possible ( Figure 2). Before the sampling, a ROI was automatically placed on the internal carotid artery to obtain an arterial input function curve, defined as the contrast concentration in the vessels feeding the tissue at each point in time during the contrast passage. The vessels, lymph nodes, and cystic areas within the salivary gland parenchyma were excluded on both DWI and DCE-PWI analyses. The ADC, I/t, AUC, and K(trans) values of the trapezius muscle were also obtained as control parameters.

Patients' Differentiations into Groups Based on Clinical Evaluation
The severity of dry mouth was clinically assessed 1 month after the end of RT by a radiotherapist with 10-years' experience using the U.S. National Cancer Institute's CTCAE v4. This score scale has three increasing levels of severity [39]. Accordingly, we Patients who were not studied with DWI and/or DCE-PWI MRI both before and after RT (n.21) Patients undergoing RT for oropharynx cancer and MRI examinations with DWI and DCE-PWI sequences both before and after the treatment (n.31) Patients who developed clinically significant xerostomia after RT and carried out MRI with DWI and DCE-PWI sequences both before and after the treatment (n.27) from the native DCE-PWI images by drawing a ROI including the largest gland section possible ( Figure 2). Before the sampling, a ROI was automatically placed on the internal carotid artery to obtain an arterial input function curve, defined as the contrast concentration in the vessels feeding the tissue at each point in time during the contrast passage. The vessels, lymph nodes, and cystic areas within the salivary gland parenchyma were excluded on both DWI and DCE-PWI analyses. The ADC, I/t, AUC, and K(trans) values of the trapezius muscle were also obtained as control parameters.  The following morphologic, DWI, and DCE-PWI features were assessed: 1. T2 signal intensity (SI) hyper-, iso-, or hypointense with respect to the muscle signal of the parotid and submandibular glands before and after RT; 2.
SI hyper-or hypointense of the parotid and submandibular glands before and after RT on DWIb800 images; 3.
Mean ADC values of the parotid and submandibular glands before and after RT (ADC pre-post) on DWI sequences; 4.
Mean AUC and K(trans) values of the parotid and submandibular glands before RT (AUCpre, K(trans)pre) and after RT (AUCpost, K(trans)post) on DCE-PWI sequences; 5.
Ratio between AUC values of parotid and submandibular glands before and after RT (AUCpost/pre); 6.
Ratio between K(trans) values of the parotid and submandibular glands before and after RT (K(trans)post/pre).
The MRI acquisition parameters are shown in Table A1 in Appendix A.

Texture Analysis
The MRI images obtained with T1 post-contrast sequences and the ADC maps before and after RT were imported into 3D slicer (www.3dslicer.org (accessed on 12 February 2022)) v10.4.2 software. The parotid and submandibular glands located on the same side of the oropharyngeal cancer, corresponding to the irradiated side, were segmented for the entirety of their volumes by a radiologist with 3-years' experience in head and neck cancer using the "segmentation wizard" extension for 3D slicer (Figure 3). Specifically, the segmentation of the gland located on the irradiated side was performed on CE-T1 sequences both before and after RT. The same method was carried out on the ADC maps thus resulting in a total of 4 different segmented volumes per gland. This process resulted in a total number of 216 volumes being segmented with 4 submandibular and 4 parotid gland volumes investigated for each of the 27 patients. Texture features were analyzed and extracted from such volumes using the extension "Pyradiomics" for 3D slicer. and after RT were imported into 3D slicer (www.3dslicer.org (accessed on 12 February 2022)) v10.4.2 software. The parotid and submandibular glands located on the same side of the oropharyngeal cancer, corresponding to the irradiated side, were segmented for the entirety of their volumes by a radiologist with 3-years' experience in head and neck cancer using the "segmentation wizard" extension for 3D slicer (Figure 3). Specifically, the segmentation of the gland located on the irradiated side was performed on CE-T1 sequences both before and after RT. The same method was carried out on the ADC maps thus resulting in a total of 4 different segmented volumes per gland. This process resulted in a total number of 216 volumes being segmented with 4 submandibular and 4 parotid gland volumes investigated for each of the 27 patients. Texture features were analyzed and extracted from such volumes using the extension "Pyradiomics" for 3D slicer.  A total of 107 radiomics features were extracted, belonging to the following categories: First Order, Shape-based (2D and 3D), Gray Level Co-occurrence Matrix, Gray Level Size Zone Matrix, Gray Level Run Length Matrix, Gray Level Dependence Matrix, and Neighboring Gray Tone Difference Matrix. The detailed explanations of the texture subclasses can be found in Table A2 (Appendix A).

Statistical Analysis
Texture analysis and functional DWI/DCE-PWI MRI were assessed with the same statistical methods.
First, a Shapiro-Wilk test was performed to determine the nature of the distribution of the data obtained from both the functional MRI and texture analysis. After separating the subjects into three groups (group 1, 2, and 3) according to their CTCAE v4 score, each parameter from the functional MRI and texture analysis was evaluated for its variation between pre-and post-RT for each group separately in order to study the correlation between RT-induced variation and the severity of xerostomia. More specifically, the parameters were tested for variations for each gland separately.
The parameters that showed normal and non-normal distributions were analyzed with the parametric t-Student test and Wilcoxon signed-rank test, respectively. A p-value < 0.01 was defined as acceptable. Once p-values were obtained using both methods (t-Student and Wilcoxon), receiver operating characteristics (ROC) curves were produced by calculating and plotting true positive and false positive rates for each statistically significant parameter within the pre-RT population to determine the area under the curve (AUC) and cut-off values to identify the subjects who would develop more or less severe xerostomia after RT (group 1 vs. group 2 vs. group 3).

Results
Eighteen patients and nine patients belonged to group 1 (mild xerostomia) and group 2 (moderate xerostomia), respectively. No group 3 patients (severe xerostomia) were found among the 27 patients enrolled. T and N stages at baseline for each patient are shown in Table 1. The radiation dose data are provided in Table A3 (Appendix A).

Morphological and Functional MRI
The SI of the submandibular and parotid glands did not change before and after RT compared to the muscle tissue on both T2 and DWI b800 images since the SI was hyperintense in all cases. In addition, the mean ADC values, K(trans) parameters, and AUC parameters did not show significant variation between pre-and post-RT. The p-value, mean value, median value, and standard deviation for the functional MRI parameters are shown in Tables 2 and 3.

Texture Analysis
No texture features showed statistically significant variation between pre-and post-RT on both the ADC maps and CE-MRI T1w sequences in group 1 (p-value > 0.01). The same results were found for group 2 on the CE-MRI T1w sequences (p-value > 0.01). On the contrary, in group 2, the ADC map values of the parotid and submandibular glands observed before RT were significantly different than those observed after RT (p-value < 0.01) ( Table 4).   After producing ROC curves on the pre-RT values, only the feature IMC 1 for the parotid glands showed an acceptable level of diagnostic accuracy (AUC = 0.727) (Figure 4).


Gray Level Non-uniformity Normalized (GLNN)-Gray Level Run Length Ma class.
After producing ROC curves on the pre-RT values, only the feature IMC 1 for parotid glands showed an acceptable level of diagnostic accuracy (AUC = 0.727) (Fig  4).

Discussion
Xerostomia is a common complication in patients receiving radiation therapy as a consequence of damage to the salivary glands. Morphological and structural changes in the irradiated glands can be non-invasively evaluated with MRI [32,40,41]. In addition, the indirect assessment of the microarchitecture of the salivary glands with texture analysis has been recently hypothesized to be a useful tool in the identification of the severity of RTinduced xerostomia [21,34,37,42]. Overall, predictive models employing texture analysis alongside imaging techniques have shown very promising results in the assessment of head and neck disease [17,18,34,43,44]. Therefore, interest towards artificial intelligence and its applications to imaging is steadily growing. In the context of the ever-growing academic importance of texture analysis [17,28,29,34], the present study was an effort to investigate its role in the assessment of salivary gland alterations and, specifically, RT-induced xerostomia.
While it is known that the development of RT-induced xerostomia correlates with the dose distribution to the salivary glands [45,46], the correlation between imaging techniques and specific alterations in the gland microarchitectural structure is still unclear [34,47]. The current study represented the first attempt to assess RT-induced xerostomia by taking advantage of different MRI techniques, including functional imaging-DCE-PWI and DWI-and texture analysis performed on CE-T1 sequences and ADC maps in both pre-and post-RT imaging.
No parameter for neither morphological (T2 sequences) nor functional (DWI and DCE-PWI) MRI yielded statistically significant results. Such findings were different from previous studies conducted by Juan et al. [33] on DCE-PWI and Zhang et al. [21,30,31] on DWI since they suggested a possible role for parameters obtained from these imaging techniques as biomarkers in the evaluation of RT-induced xerostomia. Possible explanations for the disagreement between our findings and those studies [21,30,31,33] could be that, in those last studies, a larger sample was enrolled and acid stimulation was used to assess the ADC values while also taking into account time-related parameters, such as time-to-peak ADC. In addition to that, the parameters obtained from DWI and DCE-PWI were correlated with the different amount of radiation dose emitted, whereas in our case, all patients enrolled received the same treatment.
Texture analysis techniques on CT images were used in previous studies to investigate RT-induced xerostomia [34][35][36]42], but they have not yet been employed on CE-MRI and ADC maps. In the present study, the CE-T1 MRI texture features did not show statistically significant correlation with the development of RT-induced xerostomia. On the contrary, the texture analysis carried out on the ADC maps yielded significant results. More specifically, three texture features for the submandibular glands (GLNN-First Order, IMC 2-First Order, and GLNN-Gray Level Run Length Matrix) and two texture features for the parotid glands (IMC 1 and IMC 2, both belonging to the Gray Level Run Length Matrix) showed significant differences between pre-and post-RT imaging. Among the aforementioned texture features, only IMC 1 showed acceptable levels of diagnostic accuracy (AUC = 0.727) for the development of moderate xerostomia when applied on the pre-RT population where, in fact, was found a significant decrease in its values in patients with moderate xerostomia as opposed to patients with mild xerostomia. While no feature except for IMC 1 yielded an acceptable diagnostic accuracy on its own, a more complex model taking all of them into account might better discriminate different degrees of severity of RT-induced xerostomia. In this same framework, the poor results obtained in both morphological and functional MRI would be useful as part of a wider analysis. IMC 1 is a second-order feature belonging to the gray level co-occurrence matrix group-values ranging from 0 to −∞ (values ≤ 0)-that represents a measure of the level of the textural complexity and tissue heterogeneity. In the current study, low values of IMC 1 in the pre-RT population, more precisely, values lower than the optimal cutoff of −1473 × 10 −11 , were associated with the development of moderate xerostomia (group 2). This finding seemed to be in line with other similar studies carried out on CT imaging, in which features relating to higher textural complexity and more heterogeneous distribution of grays correlated with xerostomia of greater severity [36,38]. A study by Nardone et al. [34] performed on planning CT postulated that more heterogeneous textures might be indicative of a higher salivary gland radiosensitivity. Irregular microarchitectural structure on histopathology-altered vascularization or loss of acinar cells replaced by adipose tissue-has an impact on the development of RT-induced xerostomia as stated by Teshima et al. [47]. In this context, the hypothesized higher textural complexity of the salivary glands with normal acinar tissue replaced by a variable amount of adipose tissue might explain the correlation between the severity of xerostomia and greatest values in features related to the textural heterogeneity, as suggested in previous studies on CT examinations by van Dijk et al. [36,38].
In addition, performing texture analysis on ADC maps and CE-T1 MRI sequences, as was done in our study, might prove extremely advantageous. The possibility of using MRI tools with the help of artificial intelligence to differentiate functional gland tissue from adipose tissue [17,19,48] intuitively suggested the possible benefit of these advanced techniques in the characterization of the gland radiosensitivity, albeit at the cost of a more complex standardization. It is well known that the voxel intensity on CT images relates to the intrinsic physical properties of a tissue; on the contrary, the voxel intensity on MRI acquisition techniques is highly dependent on machine-specific characteristics [44]. This makes quantitative assessments with radiomics more prone to variation based on hardwarespecific settings [44,49,50] and standardization as a whole more difficult. However, this is only partially true for quantitative functional MRI analyses, as those carried out in the present study on ADC and DCE-PWI; the hardware specifics are less impactful [44,50], which suggests other possible benefits to their implementation on texture analysis.
The main limitation of our study was the small sample size and, especially, the relatively smaller cohort of patients with moderate RT-induced xerostomia than patients with mild RT-induced xerostomia. It is reasonable to assume that all the features resulting as statistically significant in the current study may also yield acceptable results in terms of diagnostic accuracy when larger samples are selected. However, the limited number of cases enrolled (27 patients with oropharyngeal cancer) has to be connected to the originality of our study design. The assessment of both DCE-PWI and DWI in pre-and post-RT MRI and the process of texture analysis on both ADC maps and CE-T1 sequences required a wide variety of different examination and investigation techniques to be performed, thus sensibly reducing the number of eligible patients.
The authors are fully aware that a more elaborate model that combines all functional and textural parameters examined should be used to investigate such a complex phenomenon as RT-induced xerostomia. However, the necessity of a much larger sample required to design a performing classifier has made it impossible for authors to do that. This shortcoming was another limitation of the present study. The definition of a more complex classifier in the future would likely be better able to assess RT-induced xerostomia.
Finally, the results attained by our investigation were entirely related to the pre-RT cohort and, therefore, represented a cautionary glimpse into the possibilities of employing texture analysis techniques to predict salivary gland radiosensitivity before radiation therapy is administered. This remains a topic of discussion for further investigation in order to stratify patients according to the risk of developing xerostomia of different severity.

Conclusions
In our series, the texture analysis performed on the ADC maps of the parotid glands showed good accuracy in the assessment of the severity of RT-induced xerostomia in the pre-RT population (AUC = 0.727). The differentiation between mild and moderate RT-induced xerostomia was achieved with IMC1 (cut-off −1.473 × 10 −11 ) with 73% sensitivity, 75% specificity, and 75% diagnostic accuracy. Therefore, texture analysis should be considered a useful tool to estimate salivary gland radiosensitivity and help identify patients more prone to develop serious xerostomia before radiation therapy is administered.
No statistically significant parameter was found for both morphological and functional MRI (DCE-PWI and DWI) or for texture analysis on CE-T1.  Institutional Review Board Statement: The study was conducted according to the guidelines of the Declaration of Helsinki. This is a monocentric, retrospective, comparative study approved by the Ethical Review Board of the AOU Careggi (# 21800).
Informed Consent Statement: Written informed consent was obtained from all patients involved in the study.
Data Availability Statement: Not applicable.

Conflicts of Interest:
The authors declare no conflict of interest. Table A1. MRI acquisition parameters. Unenhanced scans included T1 and T2 sampling perfection with application-optimized contrasts using different flip angle evolution (SPACE) sequences with axial, coronal, and sagittal multiplanar reconstructions; axial T2 turbo spin echo; axial fat-saturated, echo-planar DWI spectral attenuated inversion recovery (SPAIR) with two b-values (b 50-800 s/mm 2 ) and ADC maps. Enhanced scans carried out after intravenous gadolinium contrast agent (gadobutrol, 1 mL/10 kg, flow 3 mL/s followed by 20 mL saline flush) consisted of an axial T1 turbo spin echo and axial T1 volumetric interpolated breath-hold examination (VIBE) Dixon.  Table A2. Texture features subclasses.

First-Order Statistics
Describes the distribution of voxel intensities within the image region defined by the mask through commonly used and basic metrics.

Gray Level Co-occurrence Matrix (GLCM)
Describes the second-order joint probability function of an image region constrained by the mask.

Gray Level Dependence Matrix (GLDM)
Quantifies gray level dependencies in an image. A gray level dependency is defined as the number of connected voxels within distance δ that are dependent on the center voxel.

Gray Level Size Zone (GLSZM)
Quantifies gray level zones in an image. A gray level zone is defined as the number of connected voxels that share the same gray level intensity.

Gray Level Run Length Matrix (GLRLM)
Quantifies gray level runs, which are defined as the length in number of pixels of consecutive pixels that have the same gray level value.

Neighboring Gray Tone Difference Matrix (NGTDM)
Quantifies the difference between a gray value and the average gray value of its neighbors within distance δ.