Texture analysis of routine T2 weighted fat-saturated images can identify head and neck paragangliomas - A pilot study.

PURPOSE: To evaluate the role of the first and second-order texture parameters obtained from T2-weighted fat-saturated DIXON images in differentiating paragangliomas from other neck masses, and to develop a statistical model to classify them.
METHOD: We retrospectively evaluated 38 paragangliomas, 18 nerve-sheath tumours and 14 other miscellaneous neck lesions obtained from an IRB approved study conducted between January 2016 and June 2019; using a composite gold standard of histopathology, cytology and DOTANOC PET CT (A total of 70 lesions in 63 patients). Fat-suppressed T2weighted-DIXON axial images were used. First and second-order texture-parameters were calculated from the original and filtered images. Feature selection using F-statistics and collinearity analysis provided 14 texture parameters for further analysis. Mann-Whitney-U test was used to compare between the groups and p-values were adjusted for multiple comparisons. ROC curve analysis was used to obtain optimal cut-offs.
RESULTS: A total of ten texture features were found to be significantly different between paragangliomas and non-paraganglioma lesions. Minimum from the histogram of grey levels was lower in paragangliomas with a cut off of ≤113.462 obtaining 62.9 % sensitivity and 77.27 % specificity in differentiating paragangliomas from non-paragangliomas. Logistic regression model was trained (n-49) using forward feature selection, which when evaluated on the validation set(n-21)- obtained an AUC of 0.855(95 %CI, 0.633 to 0.968) with a positive likelihood ratio of 4.545 (95 %CI, 1.298-15.923) in differentiating paragangliomas from non-paragangliomas.
CONCLUSION: Texture analysis of a routine imaging sequence can identify paragangliomas with high accuracy. Further development of texture analysis would enable better imaging workflow, resource utilisation and imaging cost reductions.

Introduction

Complex regional anatomy with surrounding critical vascular structures and multicompartmental spread of disease make neck masses a unique clinico-radiological challenge. Paragangliomas are one of the more frequently encountered non-epithelial neoplasms in the neck. Similar anatomical distribution (often in the carotid space) and non-specific symptomatology makes clinical differentiation of paragangliomas from other neck lesions challenging [[1], [2], [3]]. Percutaneous biopsy from paragangliomas can often mimic a wide range of benign and malignant head neck [[4], [5], [6]] pathologies. Further, the additional risks of adrenergic crisis [7,8] and bleeding from the hyper-vascular tumour makes preoperative biopsy of paragangliomas’ difficult [5,9]. Nuclear scans (In-111 pentetreotide, FDG-PET, gallium Ga-68 DOTATATE PET/CT and gallium Ga-68 DOTATOC) have a unique role in the work-up of paragangliomas. Over 40 % of paragangliomas are familial [8,10,11] and up to 85 % of familial cases may have multicentric disease [12,13]. Radioisotope studies confirm the molecular diagnosis of paragangliomas (offsetting the need for biopsy) and also evaluate disease multicentricity in the same scan, thus compensating the additional expenditure associated with these scans [14].

Semi-quantitative descriptors like lesion heterogeneity are commonly used in characterising neck lesions- for example, the salt and pepper appearance of paragangliomas and the target appearance of schwannomas. Texture analysis takes this semi-quantitative premise further [15,16]. It quantifies mathematically- the pixel to pixel heterogeneity of tumour appearance on imaging. Texture analysis can objectively predict histology and prognosticate tumours. Texture analysis based non-invasive tumour characterisation of paraganglioma from the point of care T2 weighted images could indicate a nuclear scan over further contrast administration, thus improving the imaging workflow. Given this background, we evaluate the role of the first order and second-order texture parameters obtained from T2 fat-saturated images in differentiating paragangliomas from other neck masses (including nerve sheath tumours NSTs) and develop a simple statistical model to classify them.

Materials and methods

Inclusion and exclusion criteria

This study retrospectively evaluated the fat-saturated DIXON T2 weighted images in a subset of patients with head and neck tumours, obtained during an institutional review board-approved study conducted for MRI evaluation of neck masses at our institute between January 2017 to June 2019. All patients with neck masses referred for imaging from the Otorhinolaryngology, Head & Neck Surgery department were included in this study, and informed consent was obtained. Patients with clinical or ultrasound proven thyroid or parathyroid masses, clinically evident multiple lesions, epithelial carcinomas with possible neck metastasis, and ultrasound evident cystic lesions were excluded (Fig. 1).

Fig. 1

Flow chart demonstrating the inclusion of patients in the study.

LN- lymph node.

A total of 70 lesions in 63 patients consisting of 38 paragangliomas (group-A) in 31 patients(7 patients had two lesions),18 NSTs in 18 patients (group-B) and 14 other lesions in 14 patients (group-C) were included (Table 1). Preoperative, post-operative tissue-based biopsy and DOTA-NOC PET/CT based diagnosis of paragangliomas using somatostatin receptor expression was considered as the composite gold standard.

Table 1

Pathology of the various lesions included in the study and their anatomical distribution in the supra and infrahyoid neck.

Paraganglioma		38 (54.298%)
Benign	33	28 CS/4 J F/ 1 MC
Malignant	5	5 CS
Nerve sheath tumour		18(25.71%)
Neurofibroma	2	2 CS
Schwannoma	16	13 CS/2 PPS/ SCF
Miscellaneous lesions		14(20.00%)
Metastatic Papillary carcinoma	2	1 PCS, 1 MC
Malignant Rhabdoid mesenchymal tumour	1	1 PPS
Aggressive lymphoma	1	1 CS
Chondroid lesion/ chondrosarcoma/ Osteochondroma with malignant transformation	3	1 IT/1 PPS/1 IT
Nasopharyngeal angiofibroma	2	2 MC
Synovial sarcoma	1	1 PVS
Sinonasal Glomangiopericytoma	1	1 PPS
High-grade vasoformative neoplasm/Monophasic synovial sarcoma	1	1 IF
Metastatic carcinoma -- unknown primary	2	1 PCS/PPS
Total	70

CS- carotid space, JF- Jugular Fossa, MC – multicompartmental, PPS- para-pharyngeal space, SCF- supraclavicular fossa, PCS- posterior cervical space, IT- infratemporal fossa, PVS- peri vertebral space.

MR imaging

All patients were imaged on a 3 T MRI scanner (Ingenia 3 T, Philips, The Netherlands) using a 16- channel neurovascular coil in the supine position. An axial T1-weighted DIXON sequence, T2-weighted DIXON sequence in sagittal, coronal and axial planes were obtained; subsequently, diffusion-weighted imaging, dynamic contrast-enhanced T1 perfusion imaging and post-contrast T1 weighted images in three planes were also acquired. The fat-saturated images obtained from the axial T2-weighted DIXON sequence was used for texture analysis (TR 2500−3500 ms; TE 90 ms; fast imaging mode: turbo spin echo, shot mode: multishot, TSE factor 21; slice thickness 4 mm; flip angle 90, number of signal averages - 1; field of view 190 × 190 mm; acquired matrix size 272 × 218; reconstructed matrix size- 320 × 320; reconstructed voxel size 0.59 × 0.59 mm: acquisition time ∼4 min).

Tumour segmentation

The fat-saturated T2 axial images were exported in DICOM format from the departmental PACS server and opened on an open-source image viewer- 3D Slicer 4.11.0 (https://download.slicer.org/). Two radiologists- A.G and S.M with 6 and 16 years of experience respectively in head neck imaging in consensus segmented the tumour The slice with the largest bulk of the tumour was selected, and a whole tumour 2D region of interest was drawn (Figs. 2a). Post-contrast sequences and T1 weighted images were consulted to avoid areas of haemorrhage and necrosis. The tumour periphery was not included to avoid volume averaging from surrounding tissues.

Fig. 2

a The T2 weighted DIXON fat-saturated images were opened, and two radiologists drew a 2dimensional ROI on the slice with the largest bulk of the tumour. Areas of haemorrhage, necrosis and the peripheral part of the tumour were avoided. Image normalisation and pixel resampling were done, followed by texture extraction from the original and filtered images. Laplacian of Gaussian and wavelet-based filtration was used. b A flow chart summarizing how texture features were obtained and further analysed.

Texture extraction

PyRadiomics open-source software(version, 2.0.1.) (https://pyradiomics.readthedocs.io/en/latest/index. html) [17] was used to extract texture features. A Laplacian of Gaussian (LoG) filter with sigma values of 2, 3, 4, and 5 mm; and wavelet transform using high- and low-frequency band-pass filter combinations were used for image filtration. Pixel thickness rescaling using cubic B-spline interpolation (resultant pixel size 2 × 2 × 2 mm3); along with grey-level normalisation using a ± 3σ technique and grey-level discretisation (bin–width, 5) was check once performed (Fig. 2 a and b). First-order texture parameters and second-order grey level co-occurrence matrix (GLCM) parameters were extracted from the original as well as the filtered images. A total of 520 texture features were obtained from each lesion (Table 2)

Table 2

First-order and second-order grey-level co-occurrence matrix-based texture parameters were obtained from the normalized original images and from filtered images. Wavelet-based and Laplacian of Gaussian Filters were used. A total of 520 texture parameters were obtained.

Informational Measure of Correlation (IMC) 1; Informational Measure of Correlation (IMC) 2; Inverse Difference Moment (IDM); Maximal Correlation Coefficient (MCC); Inverse Difference; Moment Normalized (IDMN); Inverse Difference (ID); Inverse Difference Normalized (IDN).

Image interpretation

The T1weighted and the T2 weighted DIXON images (including both the non-fat-saturated and fat-saturated images) were interpreted in consensus by two radiologists (SRM (4 years' experience) and SM(16 years' experience) blinded to clinical history and patient details. Non-consensus was resolved by discussion with a third radiologist. . Each tumour was considered separately; and for every tumour- location:[CS- carotid space, JF- Jugular Fossa, MC – multi-compartmental, PPS- para-pharyngeal space, SCF- supraclavicular fossa, PCS- posterior cervical space, IT- infra-temporal fossa, PVS- peri vertebral space]; slice number; laterality (for lesion localisation); intralesional flow voids on T2 and T1 weighted images; salt and pepper appearance on T2 weighted images; necrosis; target appearance (all with binary responses “yes” or “no”); T2 weighted signal intensity (“hyperintense”, “isointense”, “hypointense”); heterogeneity(“mild”, “moderate” or “marked”); typical displacement of internal and external carotids arteries (“yes”, “no” or “not applicable”) and presence of intralesional haemorrhage on T1-weighted images (“yes” or “no”) was recorded. Each lesion was classified on a Likert scale representing the confidence that the lesion is a paraganglioma with the following grades – “0- unlikely”; “1-uncertain”; “2-probable”.

Statistical analysis

The entire data was split using stratified random sampling in a 7:3 ratio. The larger data set was used for all statistical analysis while the smaller data set was set aside for validation of the cut-offs obtained. Ranking of the texture parameters versus categorical targets of paraganglioma (group-A) and non-paraganglioma lesions (group-B + C) was done using F statistics, and the image texture variables with p values <0.05 were selected for further analysis. Out of the 520 texture parameters, 145 were selected. The highly correlated texture parameters were subsequently considered redundant and were filtered out. For this, the Pearson’s correlation coefficient between each of the 145 texture features was determined. For each texture parameter, the count of correlated columns was established for a threshold of Pearson correlation coefficient r = 0.62. The column with the most correlated additional columns was retained, and all the other correlated columns were filtered out. This was repeated until only 14 texture features were retained for further analysis.

Box and whisker plots were used for data visualization (Fig. 3). The selected texture features were tested for normality using the Kolmogorov-Smirnov (K–S) test. The Mann- Whitney U test was used to compare the texture parameters between paragangliomas and NSTs(group A versus B), and between paragangliomas versus non-paraganglioma lesions (group A versus B + C). The false discovery rate was controlled using the Benjamini–Hochberg procedure with a corrected p-value of <0.05 being taken as significant.

Fig. 3

Box and whisker plot of the texture features found to be significantly different between non-paragangliomas and paragangliomas. The texture feature plotted is displayed along the longitudinal axis. The central box represents the values from the lower to upper quartile (25 to 75 percentile). The middle line represents the median. The horizontal line extends from the minimum to the maximum value, excluding outside and far out values which are displayed as separate points. The last two boxes represent the difference in texture features between paragangliomas and nerve sheath tumours. (Inverse Difference (ID); Informational Measure of Correlation (IMC) 2).

Receiver Operator Characteristics (ROC) curve (Fig. 4) was calculated for the texture features found to be significantly different between the groups using the testing data set. Bootstrapped Youden index was used to obtain the optimal cut-off with associated sensitivity and specificity from the testing set(n-49). This cut off was then applied to the validation dataset (n-21) and sensitivity, and specificity was determined. Logistic regression models were trained to differentiate a paraganglioma from a non-paraganglioma lesion. Multinomial logistic regression with the stochastic average gradient (SAG) solver was used. All categorical data were normalised using z-score (Gaussian) normalisation. A forward feature selection loop was used to identify the best subset of features for algorithm training. The loop starts with no feature selected and iteratively the feature that improves the model the most is added to the feature set. Additional parameter optimisation loop was run to identify the optimal learning rate/ step size strategy and the prior distribution of the coefficients in the resulting model (Uniform; Gauss or Laplace). The feature selection and parameter optimization loops were run on the training set alone with the testing set being held entirely separate from model training to avoid the peeking effect. A total of two models were trained – M1- using texture features and M2 – using both texture and radiologist identified features. The models were evaluated for reproducibility on the validation set (n = 21), and the diagnostic accuracies were determined (Fig. 2 b). All data handling and logistic regression training were done on KNIME Analytics Platform- 4.1.0 (https://www.knime.com/downloads) [18]. MedCalc Statistical Software (version 14.8.1, MedCalc) was used for statistical analysis and to generate the graphs.

Fig. 4

a Receiver Operating Characteristic (ROC) curve of the true positive rate (Sensitivity) plotted as a function of the false positive rate (100-Specificity) for different cut-off points for the texture features found to be significantly different between (a) paragangliomas and non-paraganglioma lesions and (b) between paragangliomas and nerve sheath tumours. Inverse Difference (ID); Informational Measure of Correlation (IMC) 2.

Results

A total of 70 lesions, including 38 paragangliomas, 18 nerve sheath tumours and 14 other neck lesions were included in the study. Median, 5th and 95th Percentiles for the parameters across the four groups – group-A, B, C and B + C were calculated and is presented in Table 3.

Table 3

Texture parameters selected for further analysis using two-step dimensional reduction were compared between the two sets - paragangliomas versus non-paraganglioma lesions and paraganglioma versus nerve sheath tumours using Mann- Whitney-U test with an adjusted p-value of<0.05 being taken as significant. (Inverse Difference (ID); Informational Measure of Correlation (IMC) 2).

+Texture features	Nerve sheath tumour N = 13 Group B	Paraganglioma N = 27 Group A	Other neck lesions N = 9 Group C	Non-paraganglioma lesions (n = 22) Group B + C	P-value(* is significant adjusted for multiple comparisons)
	Median (5th and 95th percentile)	Median (5th and 95th percentile)	Median (5th and 95th percentile)	Median (5th and 95th percentile)	Paraganglioma versus all non-paraganglioma lesion A versus B + C	Paraganglioma versus Nerve sheath tumours
original glcm Auto correlation	1283.227 (379.532 3185.353)	1453.998 (510.98 6552.004)	1461.841 (395.824 1780.55)	1372.534 (395.824 3072.483)
wavelet-LLH glcm Id	0.173 (0.111 0.233)	0.153 (0.11 0.213)	0.227 (0.154 0.356)	0.194 (0.131 0.292)	0.012*	0.427
log-sigma-5−0-mm-3D glcm Cluster Shade	61.688 (15639.997 1512.152)	589.161 (1197.333 6698.538)	722.119 (7494.751 692.117)	29.883 (7494.751 1246.099)	0.026*	0.199
wavelet-LHL firstorder Kurtosis	4.319 (2.707 7.908)	3.678 (2.065 5.268)	3.914 (3.195 16.553)	4.219 (2.979 8.812)	0.022*	0.071
wavelet-HLH firstorder Range	277.5 (121.051 455.191)	313.932 (179.891 589.79)	204.19 (140.822 311.87)	252.416 (140.822 385.198)	0.023*	0.411
wavelet-HHH firstorder Kurtosis	3.877 (2.365 5.522)	3.547 (2.69 5.028)	3.932 (3.017 7.338)	3.926 (3.017 6.871)	0.013*	0.103
original firstorder Minimum	256.252 (3.214 444.176)	102.733 (30.566 296.279)	138.495 (8.005 222.003)	195.153 (3.214 323.267)	0.027*	0.02*
wavelet-LLH glcm Correlation	0.627 (0.392 0.712)	0.606 (0.341 0.695)	0.656 (0.482 0.764)	0.651 (0.482 0.756)	0.035*	0.153
log-sigma-4−0-mm-3D glcm Imc2	0.992 (0.928 0.998)	0.994 (0.97 0.999)	0.978 (0.937 0.991)	0.986 (0.937 0.998)	0.03*	0.479
wavelet-HLL glcm Correlation	0.1 (0.034 0.328)	0.024 (0.132 0.148)	0.01 (0.077 0.128)	0.071 (0.035 0.26)	0.07*	0.002*
wavelet-HHH glcm Cluster Prominence	74353.029 (9065.688 186994.448)	74557.012 (11612.69 469793.981)	12502.425 (1886.475 76614.47)	21464.392 (2969.003 155408.687)	0.007*	0.199
wavelet-LHL firstorder 90Percentile	56.47 (36.786 106.52)	67.207 (31.778 117.895)	41.493 (19.396 101.81)	52.527 (30.718 101.81)	0.033*	0.189
wavelet-LLL glcm Joint Entropy	7.975 (6.405 9.943)	7.688 (5.093 9.304)	8.848 (7.062 9.382)	8.016 (6.516 9.382)	0.051	0.231
log-sigma-5−0-mm-3D firstorder Skewness	0.013 (0.511 0.551)	0.116 (0.77 0.849)	0.669 (1.515 0.534)	0.053 (1.088 0.543)	0.084	0.348

Paraganglioma versus non-paraganglioma lesions (Group A versus B + C)

Of the total 14 features compared; ten texture features were found to be significantly different between the two groups (adjusted p-value <0.05) (Table.3). The diagnostic accuracies of the cut-offs obtained for the texture parameters from the testing-set, along with the associated cut-offs are provided in Table 4. The minimum first-order-histogram values obtained from the normalized original images was lower in paragangliomas [113.14 ± 96.6- mean +-S.D] as compared to non-paraganglioma lesions [185.32 ± 113.62 mean+-S.D]. Wavelet HHH glcm Cluster Prominence had the most significant AUC (0.727) with a 81.48 % sensitivity and 59.09 % specificity in paraganglioma identification. First-order-Kurtosis obtained from the wavelet LHL filtered images was lower in paragangliomas (3.71 ± 1.09 mean +-S.D) as compared to non-paraganglioma lesions (5.3 ± 3.04 mean +-S.D) with 62.96 % sensitivity and 68.18 % specificity. Similarly, first-order kurtosis from wavelet HHH filtered images was lower in paragangliomas (3.59 ± 0.75mean +- S.D) versus non-paragangliomas (4.26 ± 1.18 mean+-S.D). Amongst the other texture parameters, wavelet LLH glcm Id had an AUC of 0.71(95CI, 0.563−0.831) with a 96.3 % sensitivity and 40.91 % specificity.(Fig. 3, Fig. 4a).

Table 4

ROC curves were plotted for the true positive rate as a function of the false positive rate at different cut-off points of the texture parameters found to be significantly different between paragangliomas versus non-paragangliomas, paragangliomas versus NSTs. (Inverse Difference (ID); Informational Measure of Correlation (IMC) 2).

Paraganglioma versus non-paraganglioma lesion	ROC curves were calculated on the testing data (n = 49) and Area under the curve, bootstrapped Youden Index and associated cut-offs and diagnostic metrics were determined along with 95 % confidence intervals.						The cut-offs' obtained from the validation set was tested on the testing set(n = 21), and associated sensitivity and specificity were tabulated.
Variable	Area under the ROC curve (AUC)	95 % Confidence intervalb	Associated criterion	95 % Confidence interval	Sensitivity	Specificity	Sensitivity	Specificity
log sigma 4 0 mm 3D glcm Imc2	0.682	0.533 0.807	>0.9922	>0.9889 >0.9982	62.96 (42.4–80.6)	72.73 (49.8–89.3)	54.55 (23.38–83.25)	80 (44.39–97.48)
log sigma 5 0 mm 3D glcm ClusterShade	0.687	0.539 0.812	>519.1965	>-1508.4452 >1512.1521	51.85 (31.9–71.3)	86.36 (65.1–97.1)	36.36 (10.93–69.21)	100 (69.15–100.00)
original firstorder Minimum	0.685	0.537 0.810	≤113.462	≤75.4513 ≤241.7985	62.96 (42.4–80.6)	77.27 (54.6–92.2)	54.55 (23.38–83.25)	60 (26.24–87.84)
wavelet HHH firstorder Kurtosis	0.707	0.560 0.828	≤3.7792	≤3.5473 ≤5.1336	74.07 (53.7–88.9)	63.64 (40.7–82.8)	45.45 (16.75–76.62)	50 (18.71–81.29)
wavelet HHH glcm ClusterProminence	0.727	0.581 0.845	>24704.5087	>4365.4508 >155408.6871	81.48 (61.9–93.7)	59.09 (36.4–79.3)	81.82 (48.22–97.72)	70 (34.75–93.33)
wavelet HLH firstorder Range	0.69	0.542 0.814	>311.8698	>252.8476 >455.1905	51.85 (31.9–71.3)	81.82 (59.7–94.8)	36.36 (10.93–69.21)	80 (44.39–97.48)
wavelet LHL firstorder 90Percentile	0.678	0.530 0.805	>63.6969	>50.1495 >106.5197	55.56 (35.3–74.5)	81.82 (59.7–94.8)	54.55 (23.38–83.25)	80 (44.39–97.48)
wavelet LHL firstorder Kurtosis	0.692	0.544 0.816	≤3.8132	≤2.6763 ≤5.2679	62.96 (42.4–80.6)	68.18 (45.1–86.1)	45.45 (16.75–76.62)	40 (12.16–73.76)
wavelet LLH glcm Correlation	0.677	0.528 0.803	≤0.6686	≤0.6531 ≤0.6957	88.89 (70.8−97.6)	45.45 (24.4–67.8)	81.82 (48.22–97.72)	30 (6.67–65.25)
wavelet LLH glcm Id	0.71	0.563 0.831	≤0.2128	≤0.1949 ≤0.2701	96.3 (81.0–99.9)	40.91 (20.7–63.6)	81.82 (48.22–97.72)	40 (12.16–73.76)
Paraganglioma versus nerve sheath tumour	ROC curves were calculated on the testing data (n = 40) and bootstrapped Youden Index was obtained with associated optimal cut off; sensitivity and specificity.						The cut off obtained from the validation set was tested on the testing set(n = 16), and associated sensitivity and specificity were calculated.
original first-order Minimum	0.729	0.566 0.857	≤202.5758	≤109.211 ≤241.7985	85.19 (66.3–95.8)	61.54(31.6–86.1)	81.82 (48.22–97.72)	60 (14.66–94.73)
wavelet LLH glcm Correlation	0.641	0.474 0.786	≤0.6686	≤0.6531 ≤0.6957	88.89 (70.8–97.6)	46.15(19.2–74.9)	81.82 (48.22–97.72)	20 (0.51–71.64)

Paraganglioma versus NSTs (group A versus B)

Two texture features were found to be significantly different between NSTs and paragangliomas when adjusted for multiple comparisons. Paragangliomas had a lower minimum greyscale value (113.14±96.6 mean+-S.D) compared to NSTs (213.85±130.03 mean +- S.D). A cut off of ≤202.5758 had 85.19 % sensitivity and 61.54 % specificity in identifying a paraganglioma from a NST (Figs. 3 – last two plots and 4 b). GLCM correlation a second-order texture parameter from the Wavelet LLH filtered images were lower in paragangliomas versus NSTs with an AUC of 0.641.

Radiological image analysis

36 paragangliomas (n = 38) were hyperintense on T2 weighted images [94.74 % (95CI, 82.25–99.36)] with 26 showing mild [68.42 % (95CI, 51.35–82.5)] and 11 demonstrating moderate T2heterogeneity[28.95 % (95CI, 15.42–45.9)]. Intralesional flow voids were seen in 23 masses on the T2 weighted images[60.53 % (95CI, 43.39–75.96)] and in 16 lesions on T1 weighted images[42.11 % (95CI, 26.31–59.18)]. The characteristic displacement of the internal and external carotids arteries was seen in a total of 31 tumours [81.58 % (95CI, 65.67–92.26)]. The classical salt and pepper appearance was seen in 13 cases[33.33 % (95CI, 20.63–49.02)], with intra lesional haemorrhage being present in 8 [21.05 % (95CI, 9.55–37.32)]. Necrosis was noted in 8 [21.05 % (95CI, 9.55–37.32)] giving rise to a spurious target appearance in 4 of them [10.53 % (95CI, 2.94–24.8)]. All the NSTs (n = 18) were hyperintense on the T2w sequences, with ten showing mild heterogeneity [55.56 % (95 %CI 30.76–78.47)] and 7 being moderately heterogeneous [38.89 % (95 %CI 17.3–64.25)]. 5 out of 18 NSTs showed intra lesional flow void on T2w images [27.78 % (95 %CI 9.69–53.48)]; with necrosis being seen in 9 [50 % (95 %CI 26.02–73.98)]; 10 in 18 of the NST had the characteristic target appearance [55.56 % (95 %CI 30.76–78.47)]. None of the NSTs showed salt and pepper appearance. All of the other pathologies (n = 14) were T2 hyperintense; with 7 [50 % (95CI, 23.04–76.96)] showing mild and 5 [50 % (95CI, 23.04–76.96)] moderate heterogeneity. Intra lesional flow void on T2 weighted sequences were present in 5 [35.71 % (95CI, 12.76–64.86)]. Necrosis was present in 7[50 % (95CI, 23.04–76.96)] with a consequent target appearance seen in 4 lesions [28.57 % (95CI, 8.39–58.1)]. The Chi square test showed that intra lesional flow void on T2weighted images (Pearson Chi-Square 6.174; p-value 0.046) and on T1 weighted images (Pearson Chi-Square 11.239; p-value 0.004); salt and pepper appearance on T2weighted images (Pearson Chi-Square 10.742; p-value 0.005); the presence of necrosis (Pearson Chi-Square 6.461; p-value 0.04); the classic target appearance (Pearson Chi-Square 13.04; p-value 0.001) and the typical displacement of internal and external carotids arteries (Pearson Chi-Square 24.437; p-value <0.001) was significantly different between the above three groups.

19 out of 38 paragangliomas were confidently identified as such (“2- probable”) [50 % (95CI, 33.38–66.62)], while 14 were classified as “1- uncertain” [36.84 % (95CI, 21.81–54.01)] and 5 as “0- unlikely” [13.16 % (95CI, 4.41–28.09)]. Out of the 18 NSTs; 14 were classified as being “0- unlikely” [77.78 % (95CI, 52.36−93.59)], and 4 [22.22 % (95CI, 6.41−47.64)]as being” 1- uncertain” (Table 5). Classification of a lesion as “2-probable” or “0-unlikely” had 79.17 %(95CI, 57.849%–92.868%) specificity and 100.00 %(87.656%–100.000%) sensitivity with AUC of 0.896 (96CI, 0.780 to 0.963) in characterising the lesion as paraganglioma. However, on conventional image analysis, 19 of a total 70 [27.14 % (95CI, 17.2–39.1) cases were tagged “1- uncertain”.

Table 5

The diagnostic performance of conventional imaging in paraganglioma identification.

	Confidence lesion is a Paragangliomas.	Paraganglioma versus Non- Paraganglioma lesions
		Count	Non-Paraganglioma lesions	Count	Paraganglioma	Row Total
			Count/Row Total%		Count/Row Total%
Full dataset (n = 70)	0- unlikely	28	84.85 % (95CI,68.1−94.89)	5	15.15% (95CI,5.11−31.9)	33
	1-uncertain	4	22.22 % (95CI,6.41−47.64)	14	77.78 % (95CI,52.36−93.59)	18
	2-probable	0	0% (95CI,0−17.65)	19	100 % (95CI,82.35−100)	19
Training set(n = 49)	0- unlikely	21	87.5% (95CI,67.64−97.34)	3	12.5% (95CI,2.66−32.36)	24
	1-uncertain	1	9.09% (95CI,0.23−41.28)	10	90.91 % (95CI,58.72−99.77)	11
	2-probable	0	0% (95CI,0−23.16)	14	100 % (95CI,76.84−100)	14
Testing set (n = 21)	0- unlikely	7	77.78 % (95CI,39.99−97.19)	2	22.22 % (95CI,2.81−60.01)	9
	1-uncertain	3	42.86% (95CI,9.9−81.59)	4	57.14% (95CI,18.41−90.1)	7
	2-probable	0	0% (95CI,0−52.18)	5	100 % (95CI,47.82−100)	5

Statistical modelling to identify paragangliomas from non-paraganglioma lesions

The model M1 using texture features alone obtained a sensitivity of 90.91 % (95 %CI, 58.722%–99.770%) with a specificity of 80.00 %(95 %CI, 44.390%–97.479%) in identifying a paraganglioma. It had an AUC of 0.855(95 %CI, 0.633 to 0.968) with a positive likelihood ratio of 4.545 (95 %CI, 1.298–15.923). Inclusion of radiologists identified parameters along with texture features (M2) resulted in even better AUC of 0.905(95 %CI, 0.696 to 0.988) with a sensitivity of 90.91 % (95 %CI, 58.722%–99.770%) and specificity of 90.00 % (95 %CI, 55.498%–99.747%) with a positive likelihood ratio of 9.091 (95 %CI, 1.403–58.913). While the AUC of M1, M2 and radiologist classification were not statistically different, M2 tended towards a better accuracy. In the test data (n = 21), both the model M1 and M2 accurately classified all “1-uncertain” lesions as either paraganglioma or non-paraganglioma lesion. The various parameters used in model training along with the features included in the models, model coefficients and statistics are provided in supplementary table1.

Discussion

This study demonstrated that lesion heterogeneity quantified using first and second-order texture parameters from a single routine baseline sequence can differentiate paragangliomas from non-paraganglioma neck lesions. While studies have evaluated the role of texture analysis in head-neck epithelial carcinomas and cervical lymph nodal pathologies, no related literature was found evaluating primary head-neck non-epithelial tumours. Texture analysis can differentiate squamous cell carcinomas from lymphomas [[19], [20], [21]], predict nodal spread [22], response to chemoradiotherapy [23] and p53 status [24] in head neck squamous cell carcinomas. Because of such promise, we hypothesized texture analysis could differentiate non-epithelial neck lesions as well. While advanced imaging sequences such as arterial spin labelling and dynamic contrast MRI have also found similar applications [25,26] texture analysiscan be used on routine sequences, allowing greater avaialbility. .

Imaging texture is a mathematical extension of visual cues such as heterogeneity used by radiologists [15,27]. For example, type 1 schwannomas, with high cellularity, nuclear palisading and Verocay bodies are relatively homogenous on T2 imaging. Type 2 schwannomas have loosely organised cellular structures with areas of myxomatous and cystic changes and hence are T2 heterogeneous with areas of T2 hyperintense cystic degeneration [28]. Paragangliomas, on the other hand, are highly vascular. The flow voids in addition to intratumoral bleeds give rise to the classical “salt and pepper” appearance with T2 imaging having 40 % sensitivity in diagnosing a paraganglioma [29]. Since T2 images are commonly used as a baseline image sequence, we hoped that texture analysis could objectify the visual cues from the T2 images, thus avoiding contrast administration.

First-order texture evaluates various features of the histogram obtained from the frequency of grey levels in the tumour. The first-order minimum would be the lowest greyscale value seen in the histogram of the tumour while skewness quantifies the degree of asymmetry in the distribution of grey values and how much it deviates from a normal distribution. Kurtosis similarly shows how "peaked the distribution is"; and whether the maximum grey values are concentrated towards the tail or the middle of the histogram [15,27]. Paragangliomas in our study had lower minimum greyscale histogram values as compared to NSTs and other neck lesions. This may be due to paragangliomas having intra-tumoural flow voids on T2weighted images (which appear dark) with the resultant lower minimum greyscale values. Paragangliomas also had a lower kurtosis on the wavelet filtered images as compared to both NSTs and other neck lesions. The salt and pepper appearance of paragangliomas would result in a less peaked distribution of the grey values in the histogram and hence the lower kurtosis values (Fig. 5). Meanwhile the T2 hyperintensity seen because of the myxomatous areas in NSTs, would result in a greater peaked distribution of the grey values. Other neck lesions are usually also more T2 hyperintense, thus explaining their more peaked histograms.

Fig. 5

Paragangliomas (lower panel) are highly vascular. Intra-tumoural flow voids on the T2weighted images would result in lower minimum greyscale values on the histogram as compared to NSTs(upper panel). NSTs are more heterogenous with Antoni A and B areas providing different signals on T2 imaging with multiple T2 hyperintense areas. This results in more significantly peaked curve of the histogram as compared to paragangliomas.

Second-order parameters describe the spatial distribution of the pixels in a tumour and would quantify the voxel to voxel variability of lesions based on the underlying differences in tissue architecture. While a detailed discussion of the underlying mathematics would be beyond the scope of this text – the several of the second-order texture parameters found significant in our study would reflect the differences in histological architecture, cellular distribution, necrosis and haemorrhage of paragangliomas and other non-paraganglioma lesions. The most significant of these were wavelet LLH glcm ID and wavelet HHH glcm Cluster Prominence with 0.71 and 0.727 respectively as the area under the ROC curves.

Texture features individually had moderate sensitivity and specificity in differentiating paragangliomas from non-paraganglioma lesions (Table 2). This was expected because of the overlap between the texture features obtained between these groups. We developed a logistic regression model, using texture parameters to classify a lesion as being either a paraganglioma or a non-paraganglioma, hoping a combination of all the texture features would allow higher accuracies. The logistic regression model M1 using texture features alone obtained a very high sensitivity with a high positive predictive value in identifying paragangliomas. The models obtained accuracy rates similar to that of paraganglioma identification by an expperienced radiologist. The additional advantage of the model-based approach would be in the classification of lesions where a definitive radiological diagnosis is uncertain (27.14 % in our study). Inclusion of radiologist opinion along with texture features tended towards even better accuracies –such a model could find utilisation as the second read to the interpreting radiologist.

Fine needle aspirates or biopsies from paragangliomas may be mistaken for neurofibromas, neurofibrosarcomas or malignant melanomas [[4], [5], [6]]. The risks of life-threatening bleed; haemorrhage induced fibrosis at the operative site and adrenergic crisis renders biopsy a less sought diagnostic pathway for paragangliomas. Imaging, thus, remains the primary workhorse of preoperative paraganglioma diagnosis. T2 imaging characteristics, classical neck space location and enhancement patterns on CEMRI are used in conjunction for the preoperative diagnosis of suspected paragangliomas with a high diagnostic accuracy [[30], [31], [32]]. However, a third of all paragangliomas are familial, and up to 85 % of familial paragangliomas are multicentric [10,11]. Thus in paragangliomas, a follow up nuclear scan is often performed after a contrast-enhanced MRI based diagnosis to rule out multicentric disease.

Texture analysis may find a place in the imaging workup of such lesions where the model may be used for lesion identification or as a second read (Fig. 6). Patients with neck lesions may be triaged using the statistical models to undergo a DOTANOC scan based on the T2 images directly, thus foregoing gadolinium administration. At the same time, lesions predicted to be a NST or non-paraganglioma lesion may undergo a CEMRI followed by a biopsy. For such a utilization, we believe an open source easily implementable texture analysis software would have greater utility [16].

Fig. 6

A flowchart demonstrating how texture analysis may find a role in routine clinical imaging of paragangliomas.

Limitations

The relatively small sample is a significant limitation. Because of a large number of texture parameters obtained; it is imperative to avoid multiple comparisons [16]. A two-step dimensional reduction using F-statistics and collinearity analysis was used because of the highly correlated texture features [33]. False discovery rate was further controlled using the Benjamini–Hochberg procedure and a corrected p-value of <0.05 was taken as significant. Additionally, all cut off obtained from the bootstrapped Youden index as well as the logistic regression was validated on a validation set to evaluate reproducibility.

The MRI signal in T2 images is relative, with variability between scanners and protocols affecting texture features obtained [[34], [35], [36]]. The cut-offs obtained from a single institute study on a single scanner may not be generalizable. To circumvent this problem, we normalised the images before calculating texture to try and improve reproducibility. Additionally, the pixel spacing was resampled using sitkBspline to 2,2,2, and LoG and wavelet filters were used (Fig. 1) to improve the stability of the extracted features. Despite these measures, our study should be considered as proof of concept and multi-institutional studies involving multiple scanners and imaging protocols would be required to develop a clinically implementable statistical model. Alternatively, because of the open-source nature of the texture analysis algorithm used in this paper, each institute could develop its own scanner specific cut-offs.

In conclusion, though our results are preliminary and a proof of concept; nonetheless we have demonstrated that texture analysis of baseline fat-saturated T2 weighted images might be used to differentiate a paraganglioma from other non-epithelial neck tumours. After proper validation across multiple institutes and scanner protocols, the statistical models could evolve for clinical practice, triaging patients for a nuclear scan, avoiding contrast administration or biopsy and thus improving imaging workflow.

CRediT authorship contribution statement

Adarsh Ghosh: Methodology, Software, Writing - original draft. Soumya Ranjan Malla: Investigation, Writing - original draft. Ashu Seith Bhalla: Conceptualization, Writing - review & editing, Supervision. Smita Manchanda: Conceptualization, Writing - review & editing, Supervision. Devasenathipathy Kandasamy: Conceptualization, Methodology, Writing - review & editing, Supervision. Rakesh Kumar: Resources.

Declaration of Competing Interest

The authors report no declarations of interest.