Radiomic signatures reveal multiscale intratumor heterogeneity associated with tissue tolerance and survival in re-irradiated nasopharyngeal carcinoma: a multicenter study

Background Post-radiation nasopharyngeal necrosis (PRNN) is a severe adverse event following re-radiotherapy for patients with locally recurrent nasopharyngeal carcinoma (LRNPC) and associated with decreased survival. Biological heterogeneity in recurrent tumors contributes to the different risks of PRNN. Radiomics can be used to mine high-throughput non-invasive image features to predict clinical outcomes and capture underlying biological functions. We aimed to develop a radiogenomic signature for the pre-treatment prediction of PRNN to guide re-radiotherapy in patients with LRNPC. Methods This multicenter study included 761 re-irradiated patients with LRNPC at four centers in NPC endemic area and divided them into training, internal validation, and external validation cohorts. We built a machine learning (random forest) radiomic signature based on the pre-treatment multiparametric magnetic resonance images for predicting PRNN following re-radiotherapy. We comprehensively assessed the performance of the radiomic signature. Transcriptomic sequencing and gene set enrichment analyses were conducted to identify the associated biological processes. Results The radiomic signature showed discrimination of 1-year PRNN in the training, internal validation, and external validation cohorts (area under the curve (AUC) 0.713–0.756). Stratified by a cutoff score of 0.735, patients with high-risk signature had higher incidences of PRNN than patients with low-risk signature (1-year PRNN rates 42.2–62.5% vs. 16.3–18.8%, P < 0.001). The signature significantly outperformed the clinical model (P < 0.05) and was generalizable across different centers, imaging parameters, and patient subgroups. The radiomic signature had prognostic value concerning its correlation with PRNN-related deaths (hazard ratio (HR) 3.07–6.75, P < 0.001) and all causes of deaths (HR 1.53–2.30, P < 0.01). Radiogenomics analyses revealed associations between the radiomic signature and signaling pathways involved in tissue fibrosis and vascularity. Conclusions We present a radiomic signature for the individualized risk assessment of PRNN following re-radiotherapy, which may serve as a noninvasive radio-biomarker of radiation injury-associated processes and a useful clinical tool to personalize treatment recommendations for patients with LANPC. Supplementary Information The online version contains supplementary material available at 10.1186/s12916-023-03164-3.

. Candidate clinical and radiomics variables for prediction of PRNN in training cohort..... 22 Table S7.Performance of the radiomics and clinical models in training and validation cohorts....... 23 Table S8.Generalizability of the radiomic signature across medical centers and MRI parameters. .24 Table S9.Generalizability of the radiomic signature across patient subgroups....................................25 Table S10.Prognostic significance of the radiomic signature in multivariate Cox analyses................26 Table S11.Core enrichment genes of the biological processes associated with the radiomic signature.determined by the standard threshold doses, DFI, and patient's performance status.The target volumes were delineated based on fused MRI and computed tomography simulation images at a slice thickness of 3 mm.A thermoplastic shell was applied for immobilization.The gross tumor volume (GTV) included all the contrast-enhanced lesions and lymph nodes.The clinical target volumes (CTV) were designed to encompass microscopic disease including the GTV plus a 5-mm to 10-mm margin and a smaller margin (<3 mm) allowed near critical intracranial structures or the spinal cord.The CTV also included the entire nasopharynx and the lymph node positive involved regions.No elective reirradiation of uninvolved regional lymph nodes was performed.Critical normal structures, including the brainstem, spinal cord, parotid glands, optic nerves, chiasm, lens, eyeballs, temporal lobes, temporomandibular joints, mandible, and hypophysis were contoured and set as organs at risk (OARs) during optimization.Additional 2-mm to 3-mm margins were added to the CTV and GTV to create the planning target volume (PTV) for setup variability and internal motion.Cisplatin-based induction chemotherapy (cisplatin plus 5-fluorouracil, paclitaxel or gemcitabine) and/or concurrent chemotherapy (intravenous cisplatin 80-100 mg/m2 every 3 weeks) was administrated at the physician's discretion, given that the role of chemotherapy in re-irradiated patients remains unclear 3 .

Supplementary A3. Patient characteristics
Patient demographics, tumor stage and treatment parameters were collected, including age, sex, body mass index (BMI), complications of hypertension and diabetes, DFI, recurrent tumor stage, pathological type, baseline hemoglobin level, administration of concurrent chemotherapy and re-radiotherapy planning indices including radiation dose, fractionation, and GTVrecurrence.Because of the variation in the dose-fraction size of re-radiotherapy, all re-irradiation doses were converted to equivalent doses in 2-Gy fraction values using the following equation: total dose * (dose per fraction + 10) / (2 + 10), where 10 = a/b ratio for head and neck cancer 2 .All patients were restaged according to the 8th American Joint Committee on Cancer Union for International Cancer Control staging system.Radiotherapy-induced toxicities were recorded and graded according to the toxicity criteria of the Radiation Therapy Oncology Group.

Supplementary A4. Diagnosis of PRNN and follow-up protocol
PRNN was diagnosed according to clinical symptoms, nasopharyngeal endoscopy, MRI, and pathological examinations (Figure S1) 4,5 .The typical endoscopic manifestations of PRNN included a denaturalized nasopharyngeal mucosa, necrotic substances coverage, and an exposed skull base.The MRI findings of PRNN included a discontinuous nasopharyngeal mucosa line, soft tissue defects without enhancement, and skull base osteonecrosis.For lesions that make differentiating PRNN and tumor relapse difficult, endoscopic biopsies or surgery were performed for pathological confirmation.Patients who relapsed prior to or concurrent with PRNN were excluded for diagnosis of PRNN.The diagnostic workflow of PRNN was shown in supplementary Figure S2.PRNN was diagnosed by professional radiation oncologists at the four medical centres.All patients were followed up according to the routine practice of the study centres for ≥12 months after re-radiotherapy or until death.Patients were evaluated every 3 months during the first year after the completion of treatment, every 6 months during the subsequent 2 years and every year thereafter.Nasopharyngeal endoscopy, MRI of the head and neck, chest radiography, abdominal sonography, and plasma Epstein Barr virus DNA measurement were routinely performed.

Supplementary A5. ROI segmentation
The recurrent tumor localized in the nasopharynx was selected as the main region-of-interest (ROI_nx) to predict PRNN because it is within the high dose zones of the two courses of radiotherapy and the area in which PRNN events occurred (Figure S1).The whole tumor including the lesions beyond the nasopharynx and invading the skull base and intracranial cavity was considered as a second ROI (ROI_total) to explore the potential predictive efficacy of the model.ROI segmentation was performed on three axial MRI sequences by professional radiation oncologists in the four medical centres (Sun R, Ou XM, Zhang Y, and Guan J) using the ITK-SNAP software (version 3.8.0;www.itksnap.org).Each segmentation was validated by an experienced radiologist (Lv XF) and senior radiation oncologists in the four medical centres (Chen QY, Hu CS, Yi JL, and Wu DH).We used different labels to differentiate the tumor lesion localized in the nasopharynx (label 1) from those invading the skull base and intracranial cavity (label 2, Figure S3).After 4 weeks, 50 patients in the training cohort were randomly selected and their ROIs were segmented again by the same oncologist (Sun R) and another professional radiation oncologist (Luo DH), to assess the intra-reader and inter-reader reproducibility of radiomic features.

Supplementary A6. Radiomic features extraction
The extracted features consisted of three types: first order statistics, shape, and texture (Gray Level

Co-occurrence Matrix [GLCM], Gray Level Dependence Matrix [GLDM], Gray Level Run-Length Matrix
[GLRLM], and Gray Level Size Zone Matrix [GLSZM]) 6 .After applying three-dimensional wavelet or Laplacian of Gaussian (LoG) transformations to the original images, the features were recalculated and included.
Wavelet transformation decomposes the data into high and low-frequency components, to detect the hidden temporal features and improve the signal-to-noise ratio of images.LoG is an edge enhancement filter with different sigma levels (1 mm, 3 mm, and 5 mm) indicating the coarseness of textures.Note that, all the radiomic features conform to Image Biomarker Standardization Initiative (IBSI) 8 .A total of 2106 features (702 features in each MRI series) were extracted for ROI analysis.

Supplementary A7. Feature selection and radiomic signature building
To avoid model overfitting and improve performance, feature selection was performed in the training cohort as follows: first, we used intra/inter-class correlation coefficients (ICCs) to assess the reproducibility of extracted radiomic features based on re-segmented data.Only stable features with ICCs of >0.7 were selected.
Second, each feature was weighed on its value in discriminating the 1-year PRNN from the non-PRNN cases using three feature ranking methods; these methods included univariate logistic regression, lasso regression and linear support vector machines (SVM) binary classifier.Only the 100 most important features in each series were used to construct a compact candidate feature list for further analysis.Third, a fine selection approach was applied to select the most representative features from the compact candidate feature list and generate candidate classification models.We applied a grid search strategy to determine the best number of features to be remained, which is denoted by K. We filtered the selected features recursively using a recursive feature elimination (RFE) approach based on the model's accuracy to predict PRNN on repeated five-fold cross validation and retained K features per MR series, for 3K features in total.To address the multi-dimensionality and possible redundancy of the radiomic features, we limited the correlation coefficient between features to <0.7 in RFE to retain the uncorrelated features.In this study, three widely-used modelling algorithms were used to generate models, e.g., multivariate logistic regression, linear SVM, and random forest.All three modelling algorithms were performed on the three compact candidate feature lists attained in the feature ranking step and two sets of ROIs to generate a total of 3×3×2 candidate models.Finally, the model that achieved the highest performance in the training cohort was selected.

Supplementary A8. Candidate clinical variables and model building
Candidate clinical predictor variables for PRNN included age, sex, BMI, DFI, recurrent T category, GTVrecurrence, re-irradiation dose, accumulated GTV dose, combination with chemotherapy, and baseline hemoglobin level; these factors were carefully selected according to their clinical significance and possible association with PRNN as shown in previous studies.Each predictor variable was assessed with its value in discriminating the 1-year PRNN from the non-PRNN cases and categorized into a binary variable to perform subgroup analyses (Table S6).Finally, four clinical factors that had the best discriminative ability in the training cohort (including sex, age, DFI, and GTVrecurrence, AUC >0.5), along with the important dosimetric factor of re-irradiation dose were used to construct the clinical model, using the optimal modelling algorithm among multivariate Logistic regression, SVM, and Random Forest.

Supplementary A9. Model performance
The performance of the radiomic signature was comprehensively assessed as follows: first, ROC curves were used to assess the prediction accuracy of the radiomic signature.Second, stratification analyses according to different clinical factors and MRI acquisition parameters were performed to assess the generalizability of the model.Third, a calibration curve was plotted to assess the calibration of the model.Fourth, decision curve analysis was conducted to evaluate the clinical utility of the model in guiding re-radiotherapy by quantifying the net benefits.Fifth, a nomogram was designed to predict PRNN in a continuous fashion depending on the value of the radiomic signature 7 .Sixth, subgroup analyses were performed based on the cutoff value in order to assess the prognostic value of the radiomic signature for (i) PRNN; (ii) G5-NFS; and (iii) OS in univariate Kaplan-Meier curves and multivariate Cox proportional hazard models adjusted for sex, age, DFI, GTVrecurrence, and re-irradiation dose.

Supplementary A10. Transcriptome sequencing
Tissue samples obtained within four weeks of pre-treatment MRI were retrieved for the 29 patients included in the radiogenomics study (SYSUCC-C) to perform standard transcriptome next-generation sequencing.Total RNA was extracted using a mirVana miRNA Isolation Kit (Ambion) following the manufacturer's protocol.RNA integrity was evaluated using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA).Samples with an RNA Integrity Number (RIN) of ≥7 were subjected to the subsequent analysis.The RNA sequencing libraries were constructed using TruSeq Stranded Total RNA with Ribo-Zero Gold and then sequenced on the Illumina sequencing platform (HiSeqTM 2500).The 150-bp paired-end reads from RNA-seq were mapped to the human reference genome (UCSC hg38) using STAR v020201.RSEM v1.3.0 was then used to perform gene expression quantification.

Supplementary A11. GSEA
Using the signal-to-noise measure in GSEA, we ranked 22,291 annotated genes by their association with the PRNN radiomic signature (high-risk group vs. low-risk group).The predefined gene sets for the GO based biological processes were acquired from the Molecular Signatures Database (MSigDB v7.4) platform and used to perform pre-ranked GSEA with the appropriate software (GSEA 4.0.3,https://www.gsea-msigdb.org).To gain further insight into the biological implications of the individual radiomic features selected to build the signature, we performed ssGSEA using the Gene Set Variation Analysis (GSVA) package of R. For individual patients, an enrichment score was calculated for each predefined gene set in GO-based biological processes.The association of the six radiomic features used to build the PRNN radiomic signature with the GO based biological processes were analyzed using Spearman's rank correlation coefficient.

Supplementary A12. Radiomics quality score
The radiomics quality score of the dataset was high and estimated at 66.67% (24 out of 36 points): -Image protocol quality -well-documented image protocols (for example, contrast, slice thickness, energy, etc.) and/or usage of public image protocols allow reproducibility/replicability.1/1 point -Multiple segmentations -possible actions are: segmentation by different physicians/algorithms/software, perturbing segmentations by (random) noise, segmentation at different breathing cycles.Analyse feature robustness to segmentation variabilities.1/1 point -Phantom study on all scanners -detect inter-scanner differences and vendor-dependent features.Analyse feature robustness to these sources of variability.1/1 point -Imaging at multiple time points -collect images of individuals at additional time points.Analyse feature robustness to temporal variabilities (for example, organ movement, organ expansion/ shrinkage).0/1 point -Feature reduction or adjustment for multiple testing -decreases the risk of overfitting.Overfitting is inevitable if the number of features exceeds the number of samples.Consider feature robustness when selecting features.
3/3 points -Multivariable analysis with non-radiomics features (for example, EGFR mutation) -is expected to provide a more holistic model.Permits correlating/inferencing between radiomics and non-radiomics features.1/1 point -Detect and discuss biological correlates -demonstration of phenotypic differences (possibly associated with underlying gene-protein expression patterns) deepens understanding of radiomics and biology.1/1 point -Cut-off analyses -determine risk groups by either the median, a previously published cut-off or report a continuous risk variable.Reduces the risk of reporting overly optimistic results.1/1 point -Discrimination statistics -report discrimination statistics (for example, C-statistic, ROC curve, AUC) and their statistical significance (for example, p-values, confidence intervals).One can also apply resampling method (for example, bootstrapping, cross-validation).2/2 points -Calibration statistics -report calibration statistics (for example, Calibration-in-the-large/slope, calibration plots) and their statistical significance (for example, P-values, confidence intervals).One can also apply resampling method (for example, bootstrapping, cross-validation).

Supplementary A13. Random forest model parameters and the included radiomic features
Within the training cohort, six radiomic features were selected to build the PRNN signature using the random forest model (Table S5).There are 500 trees in one forest model.The max-depth of each tree is 2. We random subsampled 75% from the total training samples to build each tree.Entropy criterion and balanced class weights were also used in forest building.

RMS
where X is the ROI with Np pixels, c a value that shifts the intensities to prevent negative values in X.This ensures that voxels with the lowest gray values contribute the least to Energy, instead of voxels with gray level intensity closest to 0.
(ii) Sphericity: shape feature, measuring the roundness of the shape of ROI relative to a sphere.It is a dimensionless measure, independent of scale and orientation.The value range is 0<sphericity≤1, where a value of 1 indicates a perfect sphere (a sphere has the smallest possible surface area for a given volume, compared to other solids).

Sphericity
Cluster Prominence (CP) of GLCM: texture feature, measuring the skewness and asymmetry of the texture.A higher value implies more asymmetry around the mean while a lower value indicates a peak near the mean value and less variation around the mean.
where P(i, j) is the (i, j)th entry in the GLCM, Ng the number of intensity levels in the image, μx the mean gray level intensity of Px, μy the mean gray level intensity of Py.

where
V is the volume of ROI in mm3 and A the surface area of ROI in mm2.(iii) Maximum 2D Diameter (Column): shape feature, the largest pairwise Euclidean distance between ROI surface voxels in the row-slice (usually the coronal) plane.(iv) Run Entropy (RE) of Gray Level Run Length Matrix (GLRLM): texture feature, measuring the uncertainty/randomness in the distribution of run lengths and gray levels.A higher value indicates more heterogeneity in the texture patterns.Ng is the number of discreet intensity values in the image, Nr the number of discreet run lengths in the image, P(i,j|θ) the run length matrix for an arbitrary direction θ, p(i,j|θ) the normalized run length matrix, defined as (, |) = P(,|) () , and ϵ an arbitrarily small positive number (≈2.2×10−16).(v) Informational Measure of Correlation 2 (IMC 2) of Gray Level Co-occurrence Matrix (GLCM): texture feature, measuring the correlation between the probability distributions of i and j (quantifying the complexity of the texture).

Figure S1 .
Figure S1.Clinical manifestations of post-radiation nasopharyngeal necrosis.(A) Location of PRNN: usually in the recurrent tumor region and within the high dose zone of first radiotherapy.The four MRI pictures showed the same patient experiencing (i) first diagnosis of NPC; (ii) cure without PRNN after first definitive radiotherapy; (iii) diagnosis of recurrence at the same site two years after first definitive radiotherapy; (iiii) occurrence of PRNN after re-radiotherapy to the recurrent tumor region.(B) Diagnosis of PRNN faced with different scenarios in clinical practice: typical type, with apparent manifestations on both MRI and endoscope; early-stage type, with denaturalized mucosa under endoscope but no tissue defects on MRI; occult type, lesion occurred in the para-pharynx region was not seen under endoscope.Therefore, MRI and endoscope are two complementary methods for diagnosis of PRNN.(C) Outcome of PRNN: occasionally curable (patient 1), mostly intractable (patient 2) and rapid progressive (patient 3).Patient 1 had early-stage mucosa necrosis after re-radiotherapy (left) and was treated with conservative therapeutics including repeated

Figure S2 .
Figure S2.The diagnostic workflow of post-radiation nasopharyngeal necrosis.A total of 761 patients with LRNPC who received curative re-radiotherapy were recruited from four medical centres and divided into training, internal validation, and external validation cohorts.All patients were followed up for clinical outcomes including PRNN, second relapse, and death after re-radiotherapy; outpatient visit, nasopharyngeal endoscopy and MRI of the head and neck were routinely performed.Overall, 355 patients were detected that suspected of PRNN during follow-up, some of whom also had a history of relapse.Among them, 94 patients who had typical PRNN imaging manifestations and no recurrent signs were diagnosed clinically as PRNN without biopsy.The other 261 patients underwent biopsy or surgery for diagnosis and treatment of PRNN; among them, 42 patients had history of relapse before PRNN, 37 patients had relapse concurrent with PRNN, and 17 patients had no evidence of PRNN, all of whom were excluded for diagnosis of PRNN.Finally, 165 patients diagnosed of PRNN after pathological confirmation.

Figure S3 .
Figure S3.Segmentation criteria of region-of-interest.Considering the complexities of structures that nasopharyngeal carcinoma habitats, different labels were used to differentiate the tumor lesion localized in the nasopharynx (A-C, label 1, green), from that invading the skull base bone and intracranial cavity (B-D, label 2, red).Segmented lesion of label 1 refers to ROI_nx; combined lesion of label 1 and 2 refers to ROI_total.For patients with synchronous nodal recurrence, additional label (A, label 3, blue) was used to delineate the lymph nodes, but it would not be used for analysis in this study.ROI, region-of-interest.

Figure S4 .
Figure S4.Overall survival of patients with and without post-radiation nasopharyngeal necrosis.(A-C) Kaplan-Meier curves showing PRNN significantly impaired overall survival in the training (A), internal validation (B), and external validation (C) cohorts.(D-F) PRNN and massive bleeding are leading causes of death for patients with PRNN in the training (D), internal validation (E), and external validation (F) cohorts.(G-I) Disease progression is the leading cause of deaths for patients without PRNN in the training (G), internal validation (H), and external validation (I) cohorts.PRNN, post-radiation nasopharyngeal necrosis.

Figure S5 .
Figure S5.Pair-wise correlations among six selected radiomic features.The correlation coefficients of Spearman test among features were all below 0.5.

Figure S6 .
Figure S6.Confusion matrix for the prediction of PRNN using radiomics model.The confusion matrix was constructed for the training, internal validation, and external validation cohorts at an optimal threshold of 0.732.

Figure S7 .
Figure S7.Nomogram predicting risk of post-radiation nasopharyngeal necrosis at different timepoints.Depending on the radiomic signature, the probability of PRNN can be indicated at a time using a nomogram.The time is defined by the color of the line (6, 12, and 24 months post re-radiotherapy).This method illustrates the value of the signature to predict PRNN in a continuous fashion as opposed to a regular Kaplan Meier that allows only discrete levels of the predictor.The PRNN probability plots were constructed with the following formula: (t)=1-e^(-H_0 (t)*e^PI), where H_0(t) is the cumulative underlying hazard function from the Cox regression and PI is the signature (predictive index).PRNN, post-radiation nasopharyngeal necrosis.

Figure S8 .
Figure S8.Prognostic significance of the radiomic signature.(A-C) Kaplan-Meier curves of grade 5 PRNN-free survival for high risk vs low risk patients stratified by the radiomic signature in the training (A), internal validation (B) and external validation (C) cohorts.(D-F) Kaplan-Meier curves of overall survival for high risk vs low risk patients stratified by the radiomic signature in the training (D), internal validation (E) and external validation (F) cohorts.HRs were estimated using univariate Cox regression analyses.PRNN, post-radiation nasopharyngeal necrosis; OS, overall survival; HR, hazard ratio; CI, confidence interval.

Table S2 . Summary of clinical outcomes of patients in training and validation cohorts.
a PRNN-related factors included massive bleeding, central nervous infection, osteomyelitis, suicide and other events owing to PRNN-related complications.bOther complications included dysphagia, radiation encephalopathy and massive bleeding without evidences of PRNN.PRNN, post-radiation nasopharyngeal necrosis; CI, confidence interval.

Table S8 . Generalizability of the radiomic signature across medical centers and MRI parameters.
MRI, Magnetic Resonance Imaging; AUC, area under the receiver operating characteristic curve; CI, confidence interval; SYSUCC, Sun Yat-sen University Cancer Center; CHCAMS, Cancer Hospital Chinese Academy of Medical Sciences; FUSCC, Fudan University Shanghai Cancer Center; NHSMU, Nanfang Hospital of Southern Medical University.