Dynamic contrast-enhanced magnetic resonance imaging-based radiomics for the prediction of progression-free survival in advanced nasopharyngeal carcinoma

To establish a multidimensional nomogram model for predicting progression-free survival (PFS) and risk stratification in patients with advanced nasopharyngeal carcinoma (NPC). This retrospective cross-sectional study included 156 patients with advanced NPC who underwent dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI). Radiomic features were extracted from the efflux rate constant (Ktrans ) and extracellular extravascular volume (Ve ) mapping derived from DCE-MRI. Least absolute shrinkage and selection operator (LASSO) Cox regression analysis was applied for feature selection. The Radscore was constructed using the selected features with their respective weights in the LASSO Cox regression analysis. A nomogram model combining the Radscore and clinical factors was built using multivariate Cox regression analysis. The C-index was used to assess the discrimination power of the Radscore and nomogram. The Kaplan–Meier method was used for survival analysis. Of the 360 radiomic features, 28 were selected (7, 6, and 15 features extracted from Ktrans , Ve, and Ktrans +Ve images, respectively). The combined Radscore k trans +Ve (C-index, 0.703, 95% confidence interval [CI]: 0.571–0.836) showed higher efficacy in predicting the prognosis of advanced NPC than Radscore k trans (C-index, 0.693; 95% CI, 0.560–0.826) and Radscore Ve (C-index, 0.614; 95% CI, 0.481–0.746) did. Multivariable Cox regression analysis revealed clinical stage, T stage, and treatment with nimotuzumab as risk factors for PFS. The nomogram established by Radscore k trans +Ve and risk factors (C-index, 0.732; 95% CI: 0.599–0.864) was better than Radscore k trans +Ve in predicting PFS in patients with advanced NPC. A lower Radscore k trans +Ve (HR 3.5584, 95% CI 2.1341–5.933), lower clinical stage (hazard ratio [HR] 1.5982, 95% CI 0.5262–4.854), lower T stage (HR 1.4365, 95% CI 0.6745–3.060), and nimotuzumab (NTZ) treatment (HR 0.7879, 95% CI 0.4899–1.267) were associated with longer PFS. Kaplan–Meier analysis showed a lower PFS in the high-risk group than in the low-risk group (p<0.0001). The nomogram based on combined pretreatment DCE-MRI radiomics features, NTZ, and clinicopathological risk factors may be considered as a noninvasive imaging marker for predicting individual PFS in patients with advanced NPC.


Introduction
Nasopharyngeal carcinoma (NPC) is an epithelial carcinoma with a distinct geographical distribution, mainly occurring in Southern China and Southeast Asia (1). More than 70% of patients are diagnosed with advanced-stage (III-IVB) NPC at presentation, owing to its nonspecific clinical symptoms and concealed anatomical location (2). The standard treatment for patients with advanced NPC is cisplatin-based concurrent chemoradiotherapy (CCRT) with or without induction chemotherapy (IC) according to National Comprehensive Cancer Network (NCCN) guidelines (3). However, locoregional recurrence and distant metastasis are the main causes of treatment failure in advanced NPC, with 5-year overall survival (OS) rates of 67-77%, and 2-year progressionfree survival (PFS) rates of 72.9% after standard treatment (4)(5)(6). Thus, accurately identifying patients at high risks for locoregional recurrence and distant metastasis before treatment can help to determine the need for more aggressive treatments. Modifying chemoradiation protocols or using antiepidermal growth factor receptor (EGFR) monoclonal antibodies and immunotherapy can reduce locoregional recurrence and distant metastasis in these patients (7,8). Therefore, predicting the high risk of poor prognosis among patients and optimizing personalized treatment strategies are critical.
The tumor-node-metastasis (TNM) staging system is widely used to predict prognosis and facilitate treatment stratification in patients with NPC. However, this system may not be sufficiently precise, as patients with the same TNM staging show different therapeutic responses and clinical outcomes (9). While some clinical and histopathological biomarkers are associated with survival in NPC patients (10,11), the clinical utility of these biomarkers is limited and unclear. Thus, an important challenge in clinical practice is defining effective and non-invasive biomarkers for prognosis to help in the selection of patients who can most benefit from specific treatments and predict the long-term therapeutic consequences.
Radiomics is a rapidly emerging field in medicine that transforms medical images into mineable high-dimensional quantitative features via a large number of automatically extracted data-characterization algorithms associated with tumor diagnosis, prognosis, and prediction of response to treatment (12,13). In NPC, radiomic analysis from multiparametric MR images has been successfully performed to predict individual PFS in advanced NPC (4,14). Moreover, the radiomics features of MR images are useful for predicting treatment response to chemoradiotherapy (15) and IC (16,17) in patients with NPC. However, most previous studies focused on conventional MR sequences. Dynamic contrast-enhanced (DCE)-magnetic resonance imaging (MRI) is an MR perfusion technique that consists of a series of rapid contrast-enhanced T1-weighted acquisitions. With proper quantitative analysis, the data provides functional information on blood flow, vascular permeability, and angiogenesis, in addition to the morphological tumor characteristics used in clinical practice (18). Quantitative DCE-MRI parameters can be easily derived for each pixel within the tumor using commercially available software and can be visualized in a separate image for each parameter (19). This makes DCE-MRI a powerful prognostic tool. Our team also confirmed that DCE parameters can distinguish hypoxiainducible factor-1a (HIF-1a), EGFR, and Ki-67 expression, which mediates prognosis in NPC tumors (20). Recent studies showed that DCE-MRI-based radiomics is more efficient than conventional MR sequences in predicting the prognosis and treatment response in malignant gliomas (21), breast cancer (22) and rectal cancer (23). To our knowledge, no published study has used DCE-MRI-based radiomics to predict individual PFS in patients with advanced NPC.
This study established a nomogram model combining radiomic features based on DCE-MRI, NTZ treatment, and clinical risk variables to predict individual PFS in patients with advanced NPC (stage III-IVB). Patients were divided into highand low-risk groups based on the nomogram model to evaluate the model efficiency. This may provide a novel tool for risk stratification and individualized therapy protocols for NPC.

Patients
The Hainan Affiliated Hospital of Hainan Medical University Institutional Review Board for Medical Ethics approved this retrospective analysis of anonymous data and waived the requirement for informed consent. Consecutive patients newly diagnosed with NPC in our hospital between April 2018 and December 2020 were included. The inclusion criteria were: (1) pathologically confirmed stage III-IV NPC; (2) application of pre-treatment DCE-MRI; (3) Karnofsky score ≥70; and (4) complete clinical and MR image data. The exclusion criteria were: (1) contraindications for MR examination; (2) poor imaging quality; (3) history of head or neck radiotherapy or chemotherapy; (4) other primary tumors; and (5) lost to follow-up.
All patients underwent MRI of the nasopharynx and neck, computed tomography (CT) scans of the chest, abdominal ultrasound or CT with enhancement, and bone scan. Patients at risk of distant metastases also underwent whole-body 18 Ffluorodeoxyglucose positron emission tomography (PET)/CT. TNM status and clinical stage were determined according to the eighth edition of the International Union Against Cancer/ American Joint Committee on Cancer (UICC/AJCC) staging system (24).
Patient baseline clinical characteristics and pathologic data were obtained from their medical records and included age, sex, Epstein-Barr virus (EBV), clinical stage, TNM stage, pathological type, HIF-1a, EGFR, Ki-67, tens of homolog deleted on chromosome ten (PTEN), and vascular endothelial growth factor (VEGF) protein expression.

Follow-up and clinical endpoint
Patients were assessed every 3 months in the first 2 years, every 6 months in years 3-5, and then annually. The follow-up period was defined as the period from therapy initiation to the last clinical visit or death. The study endpoint was PFS, which was defined as the time from the date of treatment initiation to the date of the first disease recurrence, metastasis, death from any cause, or last follow-up, whichever came first. Disease progression was confirmed by biopsy pathology and/or imaging methods such as CT, MRI, or FDG-PET/CT.
DCE-MR images were analyzed using Omni-Kinetics software (GE Healthcare, China). The extended toft linear model was selected, using the intracranial internal carotid artery as the standard arterial input function (AIF) curve. The K trans (efflux rate constant), k ep (reflux rate constant), V e (the extracellular extravascular volume), and V p (the intravascular plasma volume fraction) maps were automatically calculated.

Tumor segmentation and radiomics feature extraction
A three-dimensional (3D) region of interest (ROI) covering the whole tumor was manually drawn on the DCE-MR images using ITK-SNAP (version 3.8.0, USA, http://www.itksnap.org) by a radiologist with 5 years of head and neck radiological diagnosis experience (Wenzhu Li) and was validated by a senior radiologist with 12 years of experience (Weiyuan Huang) who was blinded to the clinical and pathological information of the patients. The ROIs were then used to extract quantitative radiomics features for analysis. Intraclass correlation coefficient (ICC) values were applied to assess the stability and robustness of radiomics feature extraction. Fifty randomly selected patients were used to test the ICC. The radiological attending physician (Wenzhu Li) and the senior physician (Weiyuan Huang) independently delineated the ROI who were blinded to the clinical and pathological information of the patients. ICC values under 0.4, between 0.4 to 0.75, and above 0.75 indicated weak, moderate, and strong agreement respectively.
The radiomic features of the lesions were extracted using AK software (Analysis Kit; GE Healthcare). A total of 360 features were extracted from K trans and V e mapping derived from DCE-MRI, including histogram features, gray-level size zone matrix (GLSZM) features, form factor features, gray-level cooccurrence matrix (GLCM) features, and run-length matrix (RLM) features. To eliminate the differences in the value scales of the extraction features, all features were normalized before selection. Each feature for all patients was normalized with Zscores, subtracting the mean value and dividing by the standard deviation. Next, dimensionality reduction was performed on the data to eliminate those with similar characteristics after calculating the Pearson correlation coefficients. When the coefficient is larger than the threshold value (currently, the default is 0.9), one of them is removed randomly. This method ensures that the similarity of features after dimensionality reduction is not high. The radiomics workflow is shown in Figure 1.

Model building and evaluation
To reduce the dimensionality and decrease redundant information, two steps were performed for feature selection. First, univariate Cox regression was performed for each feature. Features with p<0.05 were considered to be potentially associated with PFS and, thus, included in the following process. Least absolute shrinkage and selection operator (LASSO) Cox regression analysis was subsequently performed to identify the significant prognostic features, and 10-fold crossvalidation was applied for parameters perfected and over-fitting reduction. LASSO achieves this by imposing a constraint on the model parameter (l), which causes the regression coefficients of certain variables to shrink toward zero. All features with nonzero coefficients were selected in this step. A radiomics signature (Radscore) for each patient was built via a linear combination of selected features weighted by their corresponding non-zero coefficients. The Radscores of K trans images (Radscore k trans ) V e images (Radscore Ve ), and the combination of K trans and Workflow of radiomics analysis for predicting the prognosis of advanced nasopharyngeal carcinoma. The steps were: (1) Three-dimensional manual segmentation on K trans and V e images, (2) the calculation of six types of features for each patient from the defined segmentation, (3) LASSO Cox regression for feature selection and data dimension reduction, and (4) multivariate Cox regression to develop a radiomics nomogram model. Kaplan-Meier analyses were then performed to assess the prognostic value of the model. Abbreviations: VOI, volume of interest; GLSZM, gray-level size zone matrix; GLCM, gray-level co-occurrence matrix; RLM, run-length matrix; LASSO, least absolute shrinkage and selection operator; C-index, Harrell's concordance index; K-M curve, Kaplan-Meier analysis curve. V e images (Radscore k trans +Ve ) were then calculated for each patient. The clinicopathological characteristics of each patient were assessed using univariate Cox regression analysis.
To develop an optimal model, we developed a combination model (Radscore incorporating independent clinical predictors) using multivariable Cox regression analysis. Multicollinearity was assessed using variance inflation factor (VIF), and multicollinearity was considered when the VIF value exceeded 2. The VIF value was <2 indicating no multicollinearity among the predictors. Furthermore, a nomogram model was constructed based on Radscore and significant clinicopathological features to visualize the NPC prognosis. The optimal thresholds of the radiomics nomogram-defined score were identified using X-tile software (Yale University, New Haven, CT, USA). The patients were divided into high-or low-risk groups based on the score threshold. Kaplan-Meier analysis (log-rank tests) was used to analyze differences in PFS between the high-and low-risk groups. Harrell's concordance index (C-index) was used to evaluate the discriminative ability of the prognostic models (Radscore, nomogram model) based on calculating the net benefit for patients at each threshold probability. Bootstrap analyses with 1,000 resamples were used to obtain C-index statistics that were corrected for potential overfitting. The value of the C-index ranges from 0.5 to 1.0, with 0.5 indicating random chance and 1.0 indicating perfectly accurate discrimination between the predicted probability and actual outcome.

Statistical analysis
All statistical analyses were performed using R software (version 3.5.2, https://www.rproject.org). Univariate and multivariate Cox regression was performed using the "stats" package to identify independent prognostic factors. The "glmnet" package was used for LASSO Cox regression. The "survival" package was used for Kaplan-Meier analysis (log-rank tests) analyses. The "Rms" package was used to build the nomogram models. Harrell's concordance index (C-index) was used to evaluate the nomogram models. X-tile (Yale University, New Haven, CT, USA) software was used to determine the optimal cut-off values of the radiomics nomogram-defined score. Statistical significance was set at p<0.05.

Patient characteristics
A total of 217 consecutive patients with newly diagnosed stages III-IV NPC were identified in our hospital within the selected time frame. Sixty-one patients were excluded for images degraded by an artifact (n=14), incomplete immunohistochemical results (n=12), incomplete treatment (n=15), second primary tumor (n=6), and loss to follow-up (n=14). Thus, the analysis included 156 patients. Their detailed clinical information, histological characteristics, and treatment protocols are shown in Table 1. The median follow-up time was 12.4 months (range, 2-32 months) and the median PFS was 11.5 months (range, 2-29 months). During follow-up, 23 patients developed disease progression, including 6 local progression, 6 local recurrence, 3 distant metastases and 8 deaths. Feature selection and Radscore construction The ICC value of the K trans and V e maps was 0.990 (0.957, 0.994) and 0.993 (0.977, 0.996), respectively [median (lower quartile, upper quartile)]. Of the 360 radiomic features, 28 (7, 6, and 15 features extracted from K trans , V e , and K trans +V e images, respectively) were selected using the LASSO Cox regression model based on repeated 10-fold cross-validation ( Figure 2). Figure 3 shows the selected features. A radiomics signature was constructed using the selected features and their respective weights. The Radscore calculation formula consisting of these features is presented in the Supplementary Materials. The Cindex was used to evaluate the predictive accuracy (discrimination) of the Radscore, the results of which are shown in Table 2

Nomogram construction and evaluation
Univariate Cox regression analyses identified clinical stage, T-stage, and NTZ treatment as significant prognostic factors.
The VIF values of clinical stage, T stage, NTZ treatment and Radscore k trans +Ve was 1.652, 1.674, 1.130 and 1.099 respectively, indicating no multicollinearity among the predictors. Therefore, three Radscores combined with T-stage, clinical stage, and NTZ treatment were utilized to develop the prediction models. The C-indices for the models are shown in Table 2. The Radscore k trans +Ve -clinical model yielded a maximum C-index of 0.732. Therefore, we developed the Radscore k trans +Ve -clinical nomogram (Figure 4). A lower Radscore k trans +Ve (HR 3.5584, 95%CI 2.1341-5.933), lower clinical stage (HR 1.5982, 95%CI 0.5262-4.854), lower T stage (HR 1.4365, 95%CI 0.6745-3.060), and receiving NTZ treatment (HR 0.7879, 95%CI 0.4899-1.267) were associated with longer PFS ( Figure 5). Moreover, Radscore k trans +Ve had higher classification contributions compared to clinical variables in building the combined model. X-tile software identified an optimal Radscore k trans +Veclinical nomogram cutoff score of 3.1 for PFS prediction. Based on this threshold, the patients were assigned to high-risk (n=38, 24.36%) and low-risk (n=118, 75.64%) groups. Kaplan-Meier analysis showed a much lower PFS in the high-risk group than in the low-risk group (p<0.0001) ( Figure 6).

Discussion
This study developed a nomogram based on DCE-MRIbased radiomics, clinical stage, T-stage, and NTZ treatment to individually predict PFS in patients with advanced NPC. The radiomics features from the combination of K trans and V e images had better prognostic performance than those from either K trans or V e images alone. Integrating radiomics features with independent clinical factors adequately improved the predictive efficiency compared to the Radscore model. We then used the risk scores derived from the optimal model to classify patients into high-and low-risk groups, demonstrating its promise as a new biomarker for prognostic prediction in advanced NPC.
DCE-MRI allows the quantification of various vascular biomarkers that reflect the tumor microenvironment. It assesses the process of gadolinium leakage from the intravascular to extravascular compartments, which depends on multiple factors, including blood flow, vascular permeability, microvascular attenuation, and fractional volume of extracellular extravascular space (18). The K trans and V e of DCE-MRI quantitative parameters are well-established permeability parameters of tumors (25,26). K trans stands for the amount of contrast agent transferred from blood to tissue  and blood vessel permeability. V e refers to the amount of contrast agent every unit DCE-MRI allows. Previous studies reported that pre-treatment K trans and V e maps showed promise for predicting disease-specific survival and monitoring of treatment response in patients with NPC (19, 25,26). In the present study, the multivariate Cox proportional hazards model suggested that the Radscore was the only independent prognostic biomarker. Radscores were calculated only using radiomic features derived from entire tumors on the DCE parameter maps. We confirmed that DCE is an effective approach for predicting the prognosis and treatment response of patients with NPC. In addition, several studies have demonstrated that tumor heterogeneity, which is associated with aggressiveness and treatment response, can be evaluated by texture analysis on the K trans and V e maps (27)(28)(29). Their results indicate that K trans and V e are more valuable in the evaluation of tumor heterogeneity than other DCE-MRI parameters. Similar to their findings, we also observed that The x-and y-axes show the features and their coefficients, respectively. Finally, the 7, 6, and 15 features with non-zero coefficients were extracted from K trans (A), V e (B), and K trans +V e (C) images, respectively. K trans and V e showed good performance in predicting prognosis by non-invasively characterizing intra-tumor heterogeneity. Tumor heterogeneity may be associated with tumor angiogenesis, cell proliferation, necrosis, and even different tumor gene phenotypes (29). Higher tumor heterogeneity is highly associated with a poorer prognosis, which could be secondary to intrinsic aggressive biology or treatment resistance (30)(31)(32). Moreover, several studies have reported that radiomics features based on DCE-MRI may provide further insights into the heterogeneity of the tumor microvasculature and allow monitoring of pathophysiologic changes in various aspects of tumor vascular structure and functionality to help in evaluating tumor prognosis (23, 28,33). Furthermore, our results also demonstrated the significant association of DCE-MRI-based radiomics with PFS in patients with NPC. Pak et al. (34) reported similar results for glioblastoma.
Our results showed that the radiomic combined model built on K trans and V e images demonstrated better prognostic performance than the models derived from K trans or V e images alone. Radiomic models based on multiple sequences showed greater efficiency than the model based on a single sequence (17,35) as the combined model provides more morphological and functional information derived from multiple sequences rather than only a single sequence.
Subsequently, we integrated the Radscore k trans +Ve into a nomogram with clinical risk factors (clinical stage, T stage, and NTZ treatment) in patients with advanced NPC. Nomograms provide a scoring system and a visual prediction tool to help physicians rapidly evaluate individual survival posttreatment via a simple calculation in clinical practice. Yang et al.
(4)established a nomogram that included radiomics, overall stage, and other clinical factors to predict PFS in locoregionally advanced NPC, with a C-index in the validation cohort of 0.811. Consistent with their findings, we also confirmed that the Radscore k trans +Ve -clinical nomogram was highly efficient in predicting PFS in advanced NPC. The Cindex of the Radscore k trans +Ve -clinical nomogram in our study was 0.732 (95%CI: 0.599-0.864). This is the first attempt to include NTZ treatment as a risk factor in the nomogram to predict PFS in advanced NPC. In addition, our results demonstrated the association of advanced NPC treated with NTZ treatment and longer PFS.
Furthermore, we divided patients into high-and low-risk groups based on the optimal thresholds of the radiomics nomogram-defined score. The Kaplan-Meier survival curves showed a much lower PFS in the high-risk group than that in  the low-risk group, which indicated that the radiomics nomogram may contribute to the precise stratification of patients for individual therapeutic strategies in clinical practice, thereby improving the clinical outcomes of patients with advanced NPC.
This study has some limitations. First, the sample size was not large enough for validation in an independent cohort. Nevertheless, we performed bootstrap validation, which provides a stringent assessment, and the model demonstrated good accuracy for predicting PFS. Second, the mean follow-up period was relatively  Kaplan-Meier analyses were performed to estimate progression-free survival in clinical subgroups. The patients were divided into high-and low-risk groups according to the threshold value of 3.1. Kaplan-Meier analysis showed a much lower PFS in the high-risk group than that in the low-risk group (p < 0.0001).
short and a longer follow-up period was required to predict the 5year PFS rates. Third, based on the NCCN guidelines, patients were treated with radiotherapy alone or CCRT and chemotherapy using varying doses of radiotherapy. This might be a confounding factor in the evaluation of PFS. Finally, this study did not further predict treatment response to NTZ in patients with advanced NPC due to the imbalance of sample size in the prognosis grouping. Future NPC studies will develop a new model based on pretreatment radiomics features to predict the survival benefit of NTZ for individual patients with NPC.

Conclusion
The present study developed a nomogram to estimate individual PFS in patients with advanced NPC (stage III-IVb). The nomogram successfully stratified individual patients into highand low-risk groups with distinguishable prognoses, which may provide a non-invasive and effective tool for prognostic prediction and risk stratification. Moreover, the nomogram also provided additional information on the association between longer PFS in patients with advanced NPC receiving NTZ treatment.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

Ethics statement
This study was reviewed and approved by The Hainan Affiliated Hospital of Hainan Medical University Institutional Review Board for Medical Ethics approved this retrospective analysis of anonymous data and waived the requirement for informed consent. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.

Author contributions
WL, GW, FC and WH collected and analyzed the literature, drafted the figures and wrote the paper. TL, GD and QY performed the data collection. YL performed the statistical analysis. All authors contributed to the article and approved the submitted version.