Radiomic Features Associated with Lymphoma Development in the Parotid Glands of Patients with Primary Sjögren’s Syndrome

Simple Summary Patients diagnosed with primary Sjögren’s syndrome are characterized by an increased accumulation of mucosa-associated lymphoid tissue in the salivary and lacrimal glands due to chronic inflammation. Consequently, these patients present up to 40-fold higher risk of developing lymphoma, especially in the parotid gland, compared to the healthy population. Radiomics has recently proved its value in assessing tissue heterogeneity and proposing textural features that might become surrogates for biopsy. This retrospective study aimed to assess the potential value of radiomics in discovering textural analysis biomarkers associated with lymphoma development in the parotid glands of patients with primary Sjögren’s syndrome based on MR images, which might provide new directions in assessing the disease. Abstract Non-Hodgkin Lymphoma (NHL) represents a severe complication and the main cause of morbidity in patients with primary Sjögren’s syndrome (pSS). This study aimed to assess the role of textural analysis (TA) in revealing lymphoma-associated imaging parameters in the parotid gland (PG) parenchyma of patients with pSS. This retrospective study included a total of 36 patients (54.93 ± 13.34 years old; 91.6% females) diagnosed with pSS according to the American College of Rheumatology and the European League Against Rheumatism criteria (24 subjects with pSS and no lymphomatous proliferation; 12 subjects with pSS and NHL development in the PG, confirmed by the histopathological analysis). All subjects underwent MR scanning between January 2018 and October 2022. The coronal STIR PROPELLER sequence was employed to segment PG and perform TA using the MaZda5 software. A total of 65 PGs underwent segmentation and texture feature extraction (48 PGs were included in the pSS control group, and 17 PGs were included in the pSS NHL group). Following parameter reduction techniques, univariate analysis, multivariate regression, and receiver operating characteristics (ROC) analysis, the following TA parameters proved to be independently associated with NHL development in pSS: CH4S6_Sum_Variance and CV4S6_Inverse_Difference_Moment, with an area under ROC of 0.800 and 0.875, respectively. The radiomic model (resulting by combining the two previously independent TA features), presented 94.12% sensitivity and 85.42% specificity in differentiating between the two studied groups, reaching the highest area under ROC of 0.931 for the chosen cutoff value of 1.556. This study suggests the potential role of radiomics in revealing new imaging biomarkers that might serve as useful predictors for lymphoma development in patients with pSS. Further research on multicentric cohorts is warranted to confirm the obtained results and the added benefit of TA in risk stratification for patients with pSS.


Introduction
Primary Sjögren's syndrome (pSS) represents an inflammatory autoimmune disease, commonly affecting the salivary and lacrimal glands, resulting in the progressive destruction of the glandular parenchyma, which is replaced by fat and fibrous tissue. Consequently, patients with pSS suffer from exocrine dysfunction and sicca syndrome [1,2].
The hallmark of the pathological process in patients with pSS is the accumulation of mucosa-associated lymphoid tissue (MALT) due to chronic inflammation [3]. This tissue represents a substrate for non-Hodgkin lymphoma (NHL) development in this group of patients (especially the NHL-MALT histologic subtype). It is estimated that patients with pSS have an increased risk of developing lymphoma, up to 40-fold higher than in the general population [4]. Therefore, predicting the pSS outcome both at disease onset and during follow-up is of paramount importance.
Clinical and biological predictors of lymphoma in pSS have been assessed, the two major ones being persistent salivary gland swelling and cryoglobulinemic vasculitis. Composite indexes/scores have also been proposed as valuable tools in predicting lymphoma [5,6]. Several imaging features proved to be associated with lymphoma development in pSS, such as the severity of main salivary gland (MSG) parenchymal destruction based on the ultrasonographic aspect [7,8], a parotid gland (PG) stiffness value > 11.5 kPA assessed using 2D shear-wave elastography [9], and an increased diffusion restriction with low apparent diffusion coefficient [10]. However, despite the obtained results, biopsy and histopathological analysis remain mandatory for lymphoma diagnosis [11].
Recently, radiomics has emerged as a rapidly evolving post-processing imaging technique that can extract high-dimensional quantitative data from radiological images and aims to surpass the limits of the subjective observational imaging assessment [12,13]. Textural analysis (TA) features proved to be correlated to tissue heterogeneity at a cellular level [14,15] and, therefore, might represent novel biomarkers that could improve diagnostic accuracy and facilitate the decision-making process for clinicians, especially considering the current age of targeted treatment and patient-customized medicine [16,17].
In the field of salivary gland imaging, the role of radiomics was mainly assessed in the following two clinical settings: oncology and radiation-induced xerostomia, while few studies focused on inflammatory pathologies [18].
Radiomic studies with a specific focus on pSS are scarce. One study proposed a radiomic-based evaluation of the pSS scoring system using MSG ultrasound images and identified the best-performing classifier (multilayer perceptron) among the considered artificial intelligence algorithms [19]. Additionally, the efficiency of deep learning algorithms was tested for the automated PG segmentation on ultrasound images [20].
Regarding the role of radiomics in pSS using MRI, one study revealed textural features in the lacrimal gland's parenchyma that were able to distinguish between patients with pSS and healthy controls [21]. Another MRI radiomic study found imaging biomarkers that could stage disease activity in pSS. The histogram derived from the PG segmentation on the ADC map provided quantitative parameters that reflected different tissue characteristics between pSS patients and healthy volunteers [22]. Although MRI is a useful technique in the local staging of malignant PG tumors and pSS-associated lymphomas of salivary and lacrimal glands [23], to the best of our knowledge, no radiomic study regarding lymphomatous proliferation in pSS has been performed so far.
NHL represents the main complication and cause of morbidity in pSS patients, and the gold standard diagnostic method remains biopsy, which is an invasive procedure with a nonnegligible risk of sampling error [23,24]. Although pSS-associated lymphomas often present an indolent evolution, these malignancies still present the risk of dissemination to other mucosal sites or organs [25].
Therefore, the aim of this radiomic study was to identify alternative, non-invasive MRI textural features in the PG parenchyma associated with NHL development in patients with pSS that might serve as prognostic biomarkers, complement biopsy, and provide new directions in assessing the disease.

Materials and Methods
This study was performed according to the Declaration of Helsinki and received approval from the Ethical Committee of the "Iuliu Hat , ieganu" University of Medicine and Pharmacy Cluj-Napoca (registration number: 43; date: 11 February 2022). Due to the study's retrospective design, all participants' informed consent was waived.

Patients and Standard Reference
This retrospective nonrandomized study was conducted on patients with previously documented pSS who underwent head and neck MRI examination between January 2018 and October 2022 to assess MSG.
The inclusion criteria were the fulfillment of the American College of Rheumatology and the European League Against Rheumatism (ACR/EULAR) diagnostic criteria (2016) for pSS [26] and age older than 18 years. One rheumatologist with 5 years of clinical experience evaluated all patients. The Schirmer's test and the unstimulated whole salivary flow (UWSF) were assessed for all subjects, and the EULAR Sjögren's Syndrome disease activity index (ESSDAI) was computed [27]. Biological analysis (anti-Ro/La autoantibodies and rheumatoid factor (RF)) was performed for each patient.
The exclusion criteria were secondary Sjögren's syndrome (in patients with systemic lupus erythematosus, rheumatoid arthritis, progressive systemic sclerosis, or mixed connective tissue disease), sialolithiasis, previous neck radiation, history of hepatic virus C infection and patients with sicca syndrome that did not fulfill the pSS criteria.
After applying the inclusion and exclusion criteria, a cohort of 36 consecutive patients was constituted.
For the statistical analysis of the clinical and biological features, subjects were divided into two groups: one included 12 patients with pSS and NHL at the time of the MRI examination (pSS NHL group), and another control group of 24 pSS patients without lymphomatous transformation (pSS control group).
The gold standard for pSS patients with NHL was the histopathological result obtained after core needle biopsy or surgery. All samples were analyzed in the same institution. All subjects in the pSS NHL group were confirmed with NHL-MALT subtype and on MRI presented solitary solid nodules or masses identified due to the increased diffusion restriction and extremely low apparent diffusion coefficients (<0.650).
Patients in the control group were followed-up at 6 months and 1 year after the MRI exam. They did not present any clinical or imaging changes (assessed with the ultrasonography of MSG) suggesting lymphomatous transformation compared to the initial evaluation.
For the radiomic analysis, a total of 65 PGs were included (representing the training dataset). PGs were divided into two groups: 48 PGs were included in the pSS control group (both PGs of one patient underwent segmentation and feature extraction), while 17 PGs were included in the pSS NHL group (seven patients with unilateral involvement, and five patients with bilateral involvement, respectively).

Image Acquisition
All patients underwent head and neck MRI exams on a 1.5 Tesla scanner (SIGNA™ Explorer, General Electric) using an eight-channel high-resolution head coil. The MRI protocol comprised the following standard sequences: axial and coronal FSE T1-WI and FSE T2-WI, coronal STIR PROPELLER (short tau inversion recovery with periodically rotated overlapping parallel lines with enhanced reconstruction), coronal SE T1-WI with fat saturation obtained after intravenous administration of gadolinium chelate, diffusion-weighted imaging with the ADC map, and heavily three-dimensional T2-WI (MRI sialography) covering all MSG.
For the TA, the coronal STIR PROPELLER sequence was used, acquired employing the following specifications: field of view 256 × 256 mm; slice thickness 3 mm; slice gap One head and neck specialized radiologist reviewed each MRI examination on a dedicated workstation (General Electric, Advantage workstation, 4.7 edition) and identified the best artifact-free slice for assessing PG on the STIR PROPELLER sequence. After anonymization, the selected images were retrieved in a DICOM format and imported into an open-source texture analysis software, MaZda 5 (Institute of Electronics, Technical University of Lodz, Lodz, Poland) [28].
To decrease image variations in brightness and contrast that might impair the natural texture of PG, the first step in the MaZda program was to apply a grey-level normalization technique to all images. The mean (µ) and standard deviation (σ) of grey levels of voxels inside the ROIs were computed, and all outlier voxels (beyond µ ± 3σ) were consequently removed. A preprocessing wavelet filter was applied, using high and low bandpass filters.
A total of 65 PGs divided into two groups (48 PGs in the pSS control group and 17 in the pSS NHL group) underwent segmentation and feature extraction.
The segmentation process for the pSS control group implied incorporating each PG into a 2D region of interest (ROI) using a semiautomatic algorithm. Firstly, a seed was placed inside the PG parenchyma, and the ROI was automatically traced following the PG contour using gradient and geometrical coordinates. When necessary, manual corrections were further applied. An example of PG segmentation is shown in Figure 1. The segmentation was performed for both PGs of one subject. the following specifications: field of view 256 × 256 mm; slice thickness 3 mm; slice gap 0.3 mm; echo time 60 ms; repetition time 3500 ms; inversion time 1900 ms; bandwidth 437.1 Hz/pixel; acquisition time 4 min.

Image Preprocessing and Segmentation
One head and neck specialized radiologist reviewed each MRI examination on a dedicated workstation (General Electric, Advantage workstation, 4.7 edition) and identified the best artifact-free slice for assessing PG on the STIR PROPELLER sequence. After anonymization, the selected images were retrieved in a DICOM format and imported into an open-source texture analysis software, MaZda 5 (Institute of Electronics, Technical University of Lodz, Lodz, Poland) [28].
To decrease image variations in brightness and contrast that might impair the natural texture of PG, the first step in the MaZda program was to apply a grey-level normalization technique to all images. The mean (µ) and standard deviation (σ) of grey levels of voxels inside the ROIs were computed, and all outlier voxels (beyond µ ± 3σ) were consequently removed. A preprocessing wavelet filter was applied, using high and low bandpass filters.
A total of 65 PGs divided into two groups (48 PGs in the pSS control group and 17 in the pSS NHL group) underwent segmentation and feature extraction.
The segmentation process for the pSS control group implied incorporating each PG into a 2D region of interest (ROI) using a semiautomatic algorithm. Firstly, a seed was placed inside the PG parenchyma, and the ROI was automatically traced following the PG contour using gradient and geometrical coordinates. When necessary, manual corrections were further applied. An example of PG segmentation is shown in Figure 1. The segmentation was performed for both PGs of one subject.
The segmentation process for the pSS NHL group was performed first by automatically tracing a 2D ROI that comprised the focal lesion corresponding to NHL using gradient and geometrical coordinates. Then, the PG parenchyma outside the lesion was manually delineated using a nonoverlapping 2D ROI at a 2 mm distance from the first 2D ROI ( Figure 2).
The 2D ROIs represented the regions in which the radiomic features were calculated.  The segmentation process for the pSS NHL group was performed first by automatically tracing a 2D ROI that comprised the focal lesion corresponding to NHL using gradient and geometrical coordinates. Then, the PG parenchyma outside the lesion was manually delineated using a nonoverlapping 2D ROI at a 2 mm distance from the first 2D ROI ( Figure 2).
The 2D ROIs represented the regions in which the radiomic features were calculated.

Feature Extraction
From each segmented ROI, MaZda software automatically extracted a total of 275 parameters belonging to six texture classes (Absolute gradient, Histogram, Co-occurrence Matrix, Run Length Matrix, Auto-regressive Model, and Wavelet transformation). Details regarding the extracted radiomic features and the computation settings of each class are shown in Table 1.

Feature Extraction
From each segmented ROI, MaZda software automatically extracted a total of 275 parameters belonging to six texture classes (Absolute gradient, Histogram, Co-occurrence Matrix, Run Length Matrix, Auto-regressive Model, and Wavelet transformation). Details regarding the extracted radiomic features and the computation settings of each class are shown in Table 1. From the total of 275 extracted features, the MaZda program allowed the selection of the most discriminative features through several preset reduction techniques. As a first step, the probability of classification error and average correlation coefficients (POE + ACC) reduction technique was applied [29,30], and a set of 10 features was generated. This algorithm, available within the MaZda program, selects features with the highest

Feature Selection and Statistical Analysis
From the total of 275 extracted features, the MaZda program allowed the selection of the most discriminative features through several preset reduction techniques. As a first step, the probability of classification error and average correlation coefficients (POE + ACC) reduction technique was applied [29,30], and a set of 10 features was generated. This algorithm, available within the MaZda program, selects features with the highest discriminative ability while being poorly correlated, thus making them suitable for building prediction models.
To assess the stability and reproducibility of the selected TA feature set after computing the POE + ACC reduction technique, 30 PGs (randomly chosen from both studied groups) underwent re-segmentation 1 month apart from the initial procedure. The same radiologist redefined another ROI and a second round of feature extraction was performed. Then, the intraobserver reproducibility of the radiomic features was assessed using the intraclass correlation coefficient (ICC).
Radiomic parameters that presented an ICC higher than 0.850 were regarded as stable, and their corresponding absolute values from the initial segmentation were considered suitable for the subsequent statistical analysis.
A univariate analysis test (Mann-Whitney U) was further performed to assess which features were best suited to discriminate between the pSS control group and the pSS-NHL group. The statistically significant level was set at a p-value lower than 0.05. All texture parameters that showed univariate analysis results above this threshold were excluded from further processing. The receiver operating characteristic (ROC) analysis was performed, with the calculation of the area under the curve (AUC) for parameters that demonstrated statistically significant results in the univariate analysis (p < 0.05).
Optimal cutoff values were chosen using a common optimization step that maximized the Youden index. Sensitivity (Se), specificity (Sp), positive likelihood ratio (+LR), and negative likelihood ratio (−LR), with their corresponding 95% confidence intervals (CI), were computed from the same data without further adjustments.
Parameters that showed statistically significant results in the univariate and ROC analysis were included in a multiple regression using the "enter" input model. The resulting features independently associated with lymphoma development in patients with pSS were used to generate a radiomic model, computed using the regression coefficients.
An overview of the radiomic workflow used in this study is offered in Figure 3.
To assess the stability and reproducibility of the selected TA feature set after computing the POE + ACC reduction technique, 30 PGs (randomly chosen from both studied groups) underwent re-segmentation 1 month apart from the initial procedure. The same radiologist redefined another ROI and a second round of feature extraction was performed. Then, the intraobserver reproducibility of the radiomic features was assessed using the intraclass correlation coefficient (ICC).
Radiomic parameters that presented an ICC higher than 0.850 were regarded as stable, and their corresponding absolute values from the initial segmentation were considered suitable for the subsequent statistical analysis.
A univariate analysis test (Mann-Whitney U) was further performed to assess which features were best suited to discriminate between the pSS control group and the pSS-NHL group. The statistically significant level was set at a p-value lower than 0.05. All texture parameters that showed univariate analysis results above this threshold were excluded from further processing. The receiver operating characteristic (ROC) analysis was performed, with the calculation of the area under the curve (AUC) for parameters that demonstrated statistically significant results in the univariate analysis (p < 0.05). Optimal cutoff values were chosen using a common optimization step that maximized the Youden index. Sensitivity (Se), specificity (Sp), positive likelihood ratio (+LR), and negative likelihood ratio (−LR), with their corresponding 95% confidence intervals (CI), were computed from the same data without further adjustments.
Parameters that showed statistically significant results in the univariate and ROC analysis were included in a multiple regression using the "enter" input model. The resulting features independently associated with lymphoma development in patients with pSS were used to generate a radiomic model, computed using the regression coefficients.
An overview of the radiomic workflow used in this study is offered in Figure 3.  This step-by-step feature selection method was used in previous texture analysis studies [31,32], and the resulting parameters demonstrated good discriminative ability.
The statistical analysis regarding the clinical and biological features implied assessing the differences between the means or medians using the independent-samples T test or Mann-Whitney U test, as necessary. The exact Fisher test was used to evaluate the association between categorical variables.
The statistical analysis was performed using the commercially available dedicated software, MedCalc version 14.8.1 (MedCalc Software, Mariakerke, Belgium).

Results
A total of 36 patients diagnosed with pSS referred to our imaging department during the study period (mean age 54.93 ± 13.34; age range 29-83) were included in this study ( Table 2). Statistically significant risk factors for NHL development in our cohort of patients with pSS were the ESSDAI score value and the positive presence of the rheumatoid factor. Subjects in the pSS NHL group presented a higher disease activity than the pSS control group, using an ESSDAI score cutoff value of 5 (p < 0.001). The overall ESSDAI score values were higher in the pSS NHL group (p < 0.001). Disease duration did not influence the lymphomatous transformation (p > 0.05). The rheumatoid factor was present in all patients in the pSS NHL group and only 62.5% of the subjects in the pSS control group (p = 0.016).
For the textural analysis, 48 PGs from the pSS control group and 17 PGs from the pSS NHL group underwent segmentation, and a total of 275 radiomic features were extracted. Following the POE+ACC reduction technique, 10 unique texture features with the highest discriminatory values between the two studied groups were selected. Seven of the 10 previously selected texture features showed statistically significant results in the univariate analysis (Table 3). All selected parameters presented high ICC values (≥0.850). One variation of Sum Variance, Run Length NonUniformity, and Sum Average were excluded from further analysis, as the p-value exceeded 0.05.
The receiver operating characteristics (ROC) analysis was further performed. Three parameters (CH4S6SumVarnc, CV4S6InvDfMom, and Perc1) presented 88.24% sensitivity in differentiating between the two groups of patients, with specificities of 64.58%, 77.08%, and 62.50%, respectively. The highest area under the curve (AUC) was reached by CV4S6InvDfMom (0.875). The extended ROC analysis results are presented in Table 4.
The multivariate analysis showed a coefficient of determination (R 2 ) of 0.5524, an adjusted R 2 of 0.4975, a multiple correlation coefficient of 0.7433, and a residual standard deviation of 0.3140 (Table 5). Four parameters proved to be independent predictors for NHL development in patients with pSS (CH4S6SumVarnc, Perc90, Mean, CV4S6InvDfMom). A radiomic model (RM) was generated, including two independent parameters (CH4S6SumVarnc and CV4S6InvDfMom) revealed in the multivariate analysis. Two parameters (Mean, Perc90) were excluded from the model due to the high variance inflation factor (VIF), indicating multicollinearity. RM was computed using the regression coefficients (RM = −27.7065 + 0.00417CH4S6SumVarnc + 3.3534CV4S6InvDfMom).  The seven texture features that showed statistically significant results in the univariate analysis, high ICC values, and ROC analysis were included in the multivariate regression. At the cutoff value ≥ 1.556, the RM was associated with NHL development with high sensitivity and specificity (94.12% and 85.42%, respectively), presenting an AUC of 0.931. (Table 6). The areas under the ROC curve of the independent features associated with lymphoma development in pSS and the resulting radiomic model are depicted in Figure 4.

Discussion
The exocrine glands of patients with pSS are characterized by an increased content of mucosa-associated lymphoid tissue, which consequently increases the risk of develop-

Discussion
The exocrine glands of patients with pSS are characterized by an increased content of mucosa-associated lymphoid tissue, which consequently increases the risk of developing lymphoma [3].
The ability to forecast the pSS outcome on disease onset and during follow-ups is still limited despite years of research. Significant clinical and biological predictors of lymphoma in pSS proved to be persistent enlargement of the MSG (defined as lasting at least 2 months) and mixed cryoglobulinemia [3,5,33].
Currently, there are no studies to certify the role of any imaging technique as a validated predictive tool in pSS. Therefore, discovering noninvasive imaging biomarkers that might be associated with lymphoma development in pSS represents a crucial step for further clinical and applied research.
MRI represents a valuable imaging technique in pSS diagnosis, allowing both the anatomical assessment of MSG by using T1-and fat-suppressed T2-weighted images and the functional evaluation with MRI sialography, based on the spontaneously increased signal of stagnant fluids on heavily T2-weighted sequence and simultaneously signal suppression of the adjacent tissue [2]. The association of more than 5% fat areas with diminished intact parenchyma replaced by areas with increased signal in the MSG on T2-weighted images with fat saturation, together with an increased number of foci of salivary duct ectasia (≥6), reached a 96.4% sensitivity and 100% specificity for pSS diagnosis [34]. MRI sialography of PG has outstanding diagnostic performance in pSS (sensitivity and specificity ranging between 83-96% and 83-100%) [34][35][36]. MRI also proves helpful in parotid lymphoma diagnostic guidance and local staging [37,38]. According to the algorithm proposed by Jousse-Joulin et al. [11], the association of solid and cystic lesions with very low apparent diffusion coefficient values requires a biopsy, as there is high lymphoma suspicion.
Recently, radiomics proved to be a promising tool in oncological imaging, especially in diagnosing cancer, evaluating the response to therapy, or predicting prognosis [16,17].
Therefore, the present radiomic study aimed to assess the role of TA in discovering imaging biomarkers in the PG of patients with pSS that are associated with lymphoma development and could potentially become surrogates of the histopathological results.
Using fat saturation techniques, the PG parenchymal architecture was outlined by suppressing additional signals originating from interlobular fat structures [39,40]. Therefore, the STIR sequence was selected to perform PG segmentation and TA feature extraction. The STIR technique provided a more uniform fat suppression than the fat-saturated FSE technique, especially for cervical MR imaging, where the air-filled structures and complex anatomy might generate susceptibility effects and cause magnetic field inhomogeneity [41]. Moreover, the PROPELLER technique was utilized to reduce motion artifacts, improve image quality and obtain a more homogenous set of images [42].
In this study, two independent texture analysis features were found to be associated with lymphoproliferation in pSS ( Figure 5): CH4S6SumVarnc and CV4S6InvDfMom, each presenting a good AUC (>0.750).
The combination radiomic model proved to perform better than the individual parameters, reaching an AUC of 0.931.
Sum variance (CH4S6SumVarnc) and inverse difference moment (CV4S6InvDfMom) were obtained from the co-occurrence matrix. Within a given ROI, sum variance measures the deflection extent of the sum of grey-level intensity distribution from the mean greylevel value [43]. Therefore, sum variance is a parameter that reflects heterogeneity by emphasizing the deviation of neighboring grey levels from the mean in the co-occurrence matrix [44]. In our study, sum variance presented higher values in the pSS NHL group, reflecting more inhomogeneous parenchyma of the PG in comparison to the pSS control group (234.60 vs. 199.24, p < 0.001).
artifacts, improve image quality and obtain a more homogenous set of images [42].
In this study, two independent texture analysis features were found to be associated with lymphoproliferation in pSS ( Figure 5): CH4S6SumVarnc and CV4S6InvDfMom, each presenting a good AUC (>0.750).
The combination radiomic model proved to perform better than the individual parameters, reaching an AUC of 0.931. Sum variance (CH4S6SumVarnc) and inverse difference moment (CV4S6InvDfMom) were obtained from the co-occurrence matrix. Within a given ROI, sum variance measures the deflection extent of the sum of grey-level intensity distribution from the mean greylevel value [43]. Therefore, sum variance is a parameter that reflects heterogeneity by emphasizing the deviation of neighboring grey levels from the mean in the co-occurrence matrix [44]. In our study, sum variance presented higher values in the pSS NHL group, reflecting more inhomogeneous parenchyma of the PG in comparison to the pSS control group (234.60 vs. 199.24, p < 0.001).
The inverse difference is an indicator of homogeneity [45]. A wide range in levels of grey-level co-occurrences is less quantified and consequently lowers the overall feature's value. In other words, the maximum value for this feature is obtained if there is no difference in the grey levels. The inverse difference moment is conceptually similar to the inverse difference feature. However, it gives less weight to items further away from the diagonal and is also linked to homogeneity [46,47]. In the pSS NHL group, CV4S6InvDfMom proved to be lower than in the pSS control group (0.10 vs. 0.18, p < 0.0001); therefore, a large grey level variation equivalent to an increased parenchymal inhomogeneity of PG favors NHL development.
Our results show that a high PG parenchymal heterogeneity quantified by TA features is associated with NHL development. This observation agrees with other studies The inverse difference is an indicator of homogeneity [45]. A wide range in levels of grey-level co-occurrences is less quantified and consequently lowers the overall feature's value. In other words, the maximum value for this feature is obtained if there is no difference in the grey levels. The inverse difference moment is conceptually similar to the inverse difference feature. However, it gives less weight to items further away from the diagonal and is also linked to homogeneity [46,47]. In the pSS NHL group, CV4S6InvDfMom proved to be lower than in the pSS control group (0.10 vs. 0.18, p < 0.0001); therefore, a large grey level variation equivalent to an increased parenchymal inhomogeneity of PG favors NHL development.
Our results show that a high PG parenchymal heterogeneity quantified by TA features is associated with NHL development. This observation agrees with other studies that proved that severe parenchymal MSG destruction, with a consequently increased inhomogeneity assessed on ultrasonography, could be a risk factor for progression to MALT lymphoma [7][8][9].
To the best of our knowledge, this is the first study to assess the role of radiomics in depicting features associated with lymphoma development in the PG of patients with pSS, and the obtained preliminary results are promising. However, this study has several important limitations that need to be addressed.
Firstly, TA was performed in a relatively small sample of PG (n = 65) that represented the training dataset. Consequently, the small training dataset might lead to unintentional overfitting, which would hinder the generalization of the radiomic model. To address this issue, a validation dataset is warranted [48,49]. Generally, approximately 70% of the acquired dataset is used for training, and the remaining samples are used to evaluate the classifier's performance on another validation dataset [50]. We were unable to split the acquired data into a training and validation dataset, due to the limited number of observations in the pSS-NHL group (n = 17, 26% of all observations). Although pSS is one of the most common autoimmune disorders, its prevalence in the general population is still relatively low (0.06% worldwide in the general population) [51]. This fact, together with the monocentric aspect of this study, has unfortunately contributed to a limited number of patients with pSS referred to our imaging center, and an even lower number of pSS patients that developed lymphoma. Therefore, the discovered radiomic features associated with lymphoma development in the parotid glands of patients with pSS could not be tested on a separate validation set in this study. This would have significantly increased the reliability of the obtained results and counteracted any potential overfitting. More extensive multicentric studies (which would also provide external datasets) are required to collect a sufficient number of pSS subjects that could be appropriately used for training and validation groups in AI algorithms based on TA.
Moreover, important limitations in MRI radiomic studies are related to the presence of confounding factors, some related to differences in image acquisition [52][53][54][55], scanner [56], or vendor [57,58] differences. Variations in image acquisition settings such as technical parameters (matrix size, time of repetition, time of echo, signal-to-noise ratio,) or voxel size (slice thickness, pixel size) may lead to pictures of varying quality, which may impact the performance of the radiomic signatures and limits generalization [53,55].
MRI radiomic studies also present great potential in identifying predictive biomarkers in several head-neck pathology studies. However, due to the high variability in methodology, collective and accurate data assessment is limited [58].
Radiomics reporting guidelines, including Radiomics Quality Score (RQS) [59] or the Image Biomarker Standardization Initiative (IBSI) [60], proposed different approaches to conduct reproducible and generalizable radiomics studies. However, there is still a lack of consensus on how to control and reduce the effect of potential confounding factors. RQS stresses the importance of reporting exact details of the used imaging protocol, but no reliable strategies for reducing the confounding effects have been provided. Conversely, the IBSI guideline emphasizes having a pre-processing standardized algorithm for feature extraction and focuses less on limiting the confounding factors. Some studies have assessed radiomic feature robustness by using test-retest repeated scans or multiple MRI scans [61,62].
The means to control confounding factors in our study were by using a standardized MRI protocol, with fixed technical parameters for all patients with pSS that were examined in our department, and by performing image processing and computation before feature extraction. Although we strongly acknowledge the importance of feature robustness assessment that impacts model generalization [55], due to the retrospective nature of this study, additional experiments could not be performed, and therefore, we could not test the effect of confounding parameters. This represents an important limitation of the present study but represents an important objective for future prospective radiomic studies in our department.
Another limitation of the study is that the radiomic features were extracted from a 2D ROI segmentation. However, in daily clinical practice, a 3D segmentation of the PG, which has an irregular shape, might be harder to adopt, given the longer time required for segmentation and possible increased operator variability.
Moreover, differences between PGs with lymphoma and the contralateral PGs without lymphoma in the pSS NHL group could not be tested, given the low number of subjects with unilateral PG lymphomatous proliferation. Finally, since all patients in the pSS NHL group were diagnosed with MALT lymphoma, one cannot draw generalizations about other NHL subtypes.

Conclusions
This study suggests that radiomic analysis of the parotid gland's parenchyma performed on MR images has the potential to reveal new imaging biomarkers that reflect tissue heterogeneity associated with lymphoma development in patients with pSS. However, the results obtained in this study must be confirmed in larger prospective studies, using ideally multicentric cohorts, to validate the role of textural analysis in the risk stratification of patients with pSS.