Prognostic and Treatment Guiding Significance of MRI-Based Tumor Burden Features and Nodal Necrosis in Nasopharyngeal Carcinoma

We aimed to develop a nomogram integrating MRI-based tumor burden features (MTBF), nodal necrosis, and some clinical factors to forecast the distant metastasis-free survival (DMFS) of patients suffering from non-metastatic nasopharyngeal carcinoma (NPC). A total of 1640 patients treated at Sun Yat-sen University Cancer Center (Guangzhou, China) from 2011 to 2016 were enrolled, among which 1148 and 492 patients were randomized to a training cohort and an internal validation cohort, respectively. Additionally, 200 and 257 patients were enrolled in the Foshan and Dongguan validation cohorts, respectively, which served as independent external validation cohorts. The MTBF were developed from the stepwise regression of six multidimensional tumor burden variables, based on which we developed a nomogram also integrating nodal necrosis and clinical features. This model divided the patients into high- and low-risk groups by an optimal cutoff. Compared with those of patients in the low-risk group, the DMFS [hazard ratio (HR): 4.76, 95% confidence interval (CI): 3.39–6.69; p < 0.0001], and progression-free survival (PFS; HR: 4.11, 95% CI: 3.13–5.39; p < 0.0001) of patients in the high-risk group were relatively poor. Furthermore, in the training cohort, the 3-year DMFS of high-risk patients who received induction chemotherapy (ICT) combined with concurrent chemoradiotherapy (CCRT) was better than that of those who were treated with CCRT alone (p = 0.0340), whereas low-risk patients who received ICT + CCRT had a similar DMFS to those who only received CCRT. The outcomes we obtained were all verified in the three validation cohorts. The survival model can be used as a reliable prognostic tool for NPC patients and is helpful to determine patients who will benefit from ICT.


INTRODUCTION
Nasopharyngeal carcinoma (NPC) is a typical disease, of which the incidence rate is the highest in Southeast Asia (1). In the early stage, the tumor can be successfully controlled by radiotherapy alone (2), and locoregionally advanced NPC is recommended to be treated with concurrent chemotherapy (3,4). It has been reported that the 5-year overall survival (OS) has reached 85% (4,5). However, NPC is prone to recurrence and/or metastasis after certain treatments, which is the main cause of treatment failure and the major cause of mortality of NPC patients (6)(7)(8).
In addition, patients in the same clinical stage receiving similar treatments have different survival outcomes. Therefore, it is necessary to build an effective prognostic model to identify patients with a poor prognosis for whom intensive follow-up and adjuvant chemotherapy may bring about more survival benefits. Thus, many efforts have been made to identify risk factors ranging from biomarkers, such as Epstein-Barr virus (EBV) DNA and gene expression, to radiomics (9)(10)(11)(12)(13). Multiplanar magnetic resonance imaging (MRI) is one of the methods most extensively used for the precise mapping of the tumor and the accurate evaluation of NPC. In addition to precisely mapping the tumor and defining the T and N stages of NPC, attempts have been made to include information from MRI into the prognostic analysis of cancer (12,(14)(15)(16)(17).
Although tumor volume and nodal volume are commonly acknowledged prognostic tumor burden factors (18)(19)(20), NPC is an irregular tumor that shows an expansive, infiltrating, or mixed growth pattern. For the expansive growth pattern, tumor volume is highly associated with T stage, while for the infiltrating growth pattern, it is less associated. The nasopharynx is in the subjacent skull base and has a complex anatomy; invasive NPC can infiltrate the cranial nerves through structural channels such as the pterygopalatine fossa, inferior orbital fissure, ethmoid sinus, and foramen lacerum (21)(22)(23). NPC patients with small tumor volumes but large extents of invasion might have poor prognosis. In addition to tumor and nodal volume, more information is needed to precisely reflect the tumor burden. Certain prognosticrelated structural information reflecting tumor burden, such as the sectional area and vertical dimension, is not being fully utilized. At the same time, induction chemotherapy (ICT) has been shown to be efficient with low toxicity (24)(25)(26); however, not all patients benefit from ICT. A more effective prediction model is needed to identify low-risk patients to reduce overtreatment when concurrent chemoradiotherapy (CCRT) would be enough. As previous studies have revealed that nodal necrosis and nodal level are important imaging features and are independent negative prognostic factors for NPC (14,(27)(28)(29), we also collected nodal necrosis and nodal level data from MRI. Therefore, in this study, we established an MRI-based tumor burden feature (MTBF) model and developed a nomogram based on MTBF combined with nodal necrosis and some clinical factors to predict distant metastasis in NPC patients. Furthermore, we used this survival model to further explore the relationship between patients with a high risk of poor outcomes and the corresponding therapeutic schedule, which may help in making clinical decisions for individual patients suffering from NPC.

Study Design and Participants
In this retrospective study, 1640 patients with non-metastatic NPC who were treated at Sun Yat-sen University Cancer Center (Guangzhou, China) from 2011 to 2016 were enrolled. For the inclusion and exclusion criteria of the participants, see Supplementary Material p 1. A total of 1148 and 492 patients were randomized to a training cohort and an internal validation cohort, respectively, by computer-generated random numbers [random numbers from 1 to 1640 for each patient, utilizing the function "sample" in R project (version 3.5.2) with a pre-defined seed "250, " and the ratio of training cohort to internal validation cohort size is 7 to 3]. To further validate the generalizability of the model, we enrolled 200 and 257 patients into two external validation cohorts, the Foshan and Dongguan validation cohorts, respectively. The patients in these two external validation cohorts were enrolled following the same criteria. The tumor stage of the enrolled patients was determined in accordance with the 8th edition of the American Joint Committee on Cancer staging manual. This study was carried out with approval of the ethical committees of the Chinese Clinical Trial Registry (ChiECRCT20190127). We have uploaded all crucial research data to the Research Data Deposit public platform (RDDA2020001382).

Procedures
We utilized a Microsoft Excel form to collect sociodemographic information (including age, sex, native place, weight, BMI, and height) and baseline clinical data [EBV DNA before treatment, EA-IgA, VCA-IgA, LYMPH, WBC, RBC, PLT, HGB, total protein (TP), ALB, lactate dehydrogenase (LDH), C-reactive protein (CRP), ABO blood type, RH blood type, T stage, N stage, clinical stage, and therapeutic regimen]. All variables above were categorized based on routine cutoff points in clinical findings and applications. Nodal necrosis (positive vs negative) and nodal level (above the caudal edge of cricoid cartilage vs lower) were also recorded as categorical variables. Tumor burden variables (volume, maximum cross-sectional area, and vertical dimension of the primary tumor and regional lymph nodes, of which the abbreviations are Tv and Lv, Ta and La, and Td and Ld, respectively) were also recorded in the Excel form as continuous variables.
Unenhanced and enhanced head and neck MRI of NPC participants were assessed, and the execution details are presented in Supplementary Material p 1. Based on enhanced T1-weighted imaging (T1WI), three specialists in MRI contoured the margin of the primary tumor and regional lymph nodes by utilizing Medical Imaging Interaction Toolkit (MITK) software (version MITK-2016.11.0; Supplementary Figure 1). Any disagreements were resolved by a consensus. Detailed information on the calculation of Tv, Lv, Ta, La, Td, and Ld is presented in Supplementary Materials p 1-2

Treatment Methods
All patients underwent radiotherapy. Patients from our center and those in the Foshan validation cohort received intensity-modulated radiation therapy (IMRT). A total of 2.3% (6 of 257) of the patients in the Dongguan validation cohort received two-dimensional radiotherapy (2D-RT), and 97.7% (251 of 257) of the patients received IMRT. The cumulative radiation doses were 66 Gy or greater in 30-35 fractions, and treatment was delivered once daily, over 5 fractions per week. All patients received platinum-based chemotherapy, including concurrent chemotherapy and ICT. ICT consisted of cisplatin with 5-fluorouracil, taxanes, or both triweekly for two or three cycles. Concurrent chemotherapy consisting of cisplatin was administered weekly or triweekly during radiotherapy.

Follow-Up and Study Endpoints
Follow-up surveys were conducted every 3 months over the first 2 years after radical therapy, once semiannually in the third to fifth years, and once yearly thereafter. The follow-up involved the determination of plasma EBV DNA concentration and indirect or direct nasopharyngoscopy, X-ray/plain and contrast-enhanced CT imaging of the chest, sonography/plain and enhanced CT of the abdomen, and plain and enhanced MRI of the head and neck. In this study, the primary endpoint was distant metastasisfree survival (DMFS), defined as the period from the first day of treatment (radiotherapy or chemotherapy) to the first occurrence of distant metastasis. The second endpoint was progressionfree survival (PFS), defined as the period of time from the first diagnosis to locoregional failure, distant failure, or death from any cause, whichever occurred first.

Statistical Analysis
We first carried out univariate analysis on the six tumor burden variables (Tv, Lv, Ta, La, Td, and Ld), in which the p value adopted for excluding variables with least significance was 0.1 (p > 0.1). Then, we conducted multivariate Cox regression analysis with a stepwise step to screen out variables that can be used to establish a prognostic model with the coefficients weighted by the Cox model in the training cohort. Stepwise regression introduces the independent variable one by one on the condition that the independent variable is significant after the F test, which is commonly used to eliminate multicollinearity and select the "optimal" regression equation. In the Guangzhou training cohort, the optimal threshold was determined by using X-tile version 3.6.1, a software developed by Yale University, which allows us to divide patients into high-and low-MTBF groups. The thresholds were determined as the values that can generate the maximum chi-square values in the Mantel-Cox test (30).
To further estimate the effect of MTBF (high vs low), nodal necrosis, and nodal level, we also developed three clinical nomograms after univariate analysis. Considering the overlapping information between N stage and nodal level, we included nodal level into the backward multivariate Cox regression with other clinical factors to generate nomogram A. According to Akaike's information criterion (AIC), sex, age, N stage, clinical stage, and LDH were selected to generate nomogram A to forecast DMFS in the training cohort. Since relevant results show that plasma EBV DNA is a prospective biomarker in NPC, we developed nomogram B by adding plasma EBV DNA. After that, we developed nomogram C based on nomogram B by integrating MTBF and nodal necrosis.
Then, calibration curves and the concordance index (C-index, proposed by Harrell) were calculated to assess and compare the prediction performance of the nomograms. To further explore the sensitivity and specificity of the prognostic model, we performed time-dependent receiver operating characteristic (ROC) analysis and calculated the area under the curve (AUC). To explore the association of nomogram C with DMFS, we calculated the total risk score for NPC patients based on nomogram C, and the patients were separated into high-risk and low-risk groups according to a cutoff value. Risk stratification was evaluated using Kaplan-Meier analysis, and the survival curves of the high-risk and low-risk groups were compared using the log-rank test. Differences in PFS between the two groups were also assessed.
The univariate and multivariate analyses of DMFS were carried out using the survival package in R (available online). With the stopping rule of AIC, the likelihood ratio test was carried out to apply backward stepwise selection. Then, the nomograms were developed by using the rms package in R (available online), including the coefficients in the multivariable Cox regression model. In this study, R version 3.5.2 and SPSS version 22.0 were adopted for statistical analyses, and a twosided p value < 0.05 was adopted to indicate differences with statistical significance.

Patient Data and Follow-Up
In total, 2097 patients from three Chinese academic institutions were enrolled. The study flow diagram is shown in Figure 1. Among 1640 patients treated at Sun Yat-sen University Cancer Center, 1148 and 492 patients were randomized to a training cohort and an internal validation cohort, respectively. Additionally, we enrolled 200 and 257 patients into two independent external validation cohorts, the Foshan and Dongguan validation cohorts, respectively. The median followup of the combined Guangzhou cohort was 55.9 months (IQR 41.4-68.7 months), that of the Dongguan cohort was 49.2 months (IQR 37.4-60.7 months), and that of the Foshan cohort was 42.4 months (IQR 37.3-46.0 months). Up to the last follow-up, 205 (12.5%) patients in the combined Guangzhou cohort, 31 (12.1%) in the Dongguan cohort, and 19 (9.5%) in the Foshan cohort developed distant metastasis. For more information on the patients' demographic information and baseline clinical characteristics, see Table 1.

Construction of the MTBF Model
In the Guangzhou training set, the first multivariate Cox regression analysis with a stepwise step showed that five of the six tumor burden variables, namely, Lv, Ta, La, Td, and Ld, were related to DMFS. In addition, a formula was generated based on the 5 tumor burden features weighted by their respective regression coefficients (Supplementary Table 1) to calculate the tumor burden score for these patients. The formula is as follows: We generated an optimal threshold (3.2) via X-tile plots to assign the patients in the training cohort into high-and low-MTBF groups (Supplementary Figure 2). In this section, 882 (76.8%), and 266 (23.2%) patients in the training cohort, 361 (73.4%), and 131 (26.6%) patients in the Guangzhou internal validation cohort, 204 (79.4%), and 53 (20.6%) patients in the Dongguan external validation cohort, and 162 (81.0%) and 38 (19.0%) patients in the Foshan external validation cohort were divided into the high-and low-MTBF groups, respectively. The result of univariate analysis of MTBF is listed in Table 2.

Development and Validation of Nomograms to Predict Survival
All variables were assessed primarily through univariate analysis for DMFS. The predictors included MTBF (high vs low), nodal necrosis (positive vs negative), nodal level (above the caudal edge of cricoid cartilage vs lower), sex, age, N stage, clinical stage, LDH, and plasma EBV DNA (Supplementary Table 2 Table 3 and Supplementary Figure 3A). Because relevant results showed that EBV DNA is a prospective prognostic biomarker of NPC, we used sex, age, N stage, clinical stage, LDH, and EBV DNA to develop nomogram B (C-index 0.702, 95% CI: 0.658-0.746; Supplementary Table 3 and Supplementary Figure 3B). Then, we used the six clinical factors, MTBF, and nodal necrosis to develop nomogram C (C-index 0.741, 95% CI: 0.702-0.780; Figure 2A). The calibration curves for 3-year DMFS also performed well in the Guangzhou internal validation set (C-index 0.738, 95% CI: 0.676-0.780), Dongguan external validation set (0.747, 0.678-0.816), and Foshan validation set (0.757, 0.657-0.857; Figures 2B-E). Then, we compared the sensitivity and specificity of the three nomograms through ROC analysis and found that nomogram C had a higher AUC value than the other two nomograms (Figure 3 and Supplementary Table 4). The ROC curve of nomogram C vs nomogram B or A reached statistical significance in the Guangzhou data set, but that of nomogram C vs B was not statistically significant in the two external validation sets.

Nomogram C and Treatment Direction
In the Guangzhou training cohort, 486 and 439 patients in the low-risk group and 87 and 136 patients in the high-risk group received CCRT and ICT + CCRT, respectively. In the high-risk group, the DMFS of patients treated with ICT + CCRT was longer than that of those treated with CCRT alone (HR 0.60, 95% CI: 0.37-0.97, p = 0.0340; Figure 5A). However, the treatment effects of ICT + CCRT and CCRT were similar in the low-risk HGB, hemoglobin; LDH, lactate dehydrogenase; EBV DNA, Epstein-Barr virus DNA; Tv, volume of the primary tumor; Lv, volume of the regional lymph nodes; Ta, maximum cross-sectional area of the primary tumor; La, maximum cross-sectional area of the regional lymph nodes; Td, vertical dimension of the primary tumor; and Ld, vertical dimension of the reginal lymph nodes. group ( Figure 5D). This finding was validated in the internal and combined external validation sets (Figures 5B,C,E,F). Compared with patients not receiving ICT, among all the cohorts, treatment with ICT was also related to an improvement of PFS in patients with high-risk scores but not in those with low-risk scores (Supplementary Figure 4).

DISCUSSION
In this retrospective multicenter cohort study, we extracted six multidimensional tumor burden-related variables, including the volume, section area, and vertical dimension of the primary tumor and regional lymph nodes, based on which we developed the MTBF model. Compared to nomogram A (based on routine clinical features) and nomogram B (nomogram A + EBV DNA), nomogram C incorporating MTBF and nodal necrosis to nomogram B had the highest C-index (0.741, 95% CI: 0.702-0.780) and AUC value (0.761, 95% CI: 0.719-0.802). Compared with those of patients in the low-risk group stratified by nomogram C, the DMFS (HR: 4.76, 95% CI: 3.39-6.69; p < 0.0001), and PFS (HR: 4.11, 95% CI: 3.13-5.39; p < 0.0001) of patients in the high-risk group were relatively poor. Moreover, compared with patients with low risk, treatment with ICT + CCRT tended to have better effects among those with high risk. Tumor load heterogeneity in the same T and N classification could cause much difficulty for predicting prognosis. Though the volume of the primary tumor is now a commonly acknowledged tumor burden factor for prognosis prediction (18)(19)(20), the mining of other information reflecting tumor burden remains insufficient. Our findings, which extracted tumor burden profiles such as the volume, section area, and vertical dimension of the primary tumor and regional lymph nodes based on pretreatment MRI, showed that MTBF is a promising risk factor (HR 4.52, 95% CI: 3.21-6.37; p < 0.0001; Supplementary Table 2). Intriguingly, Tv was removed from the construction of the MTBF model during the stepwise regression, whereas the other five tumor burden variables were retained. These results also reflect that the maximum cross-sectional area and vertical dimension could complement the effect of volume on prognosis, which offers more detailed prognostic information. In addition, compared with nomogram A (including sex, age, N stage, clinical stage, and LDH) and nomogram B (including the variables in nomogram A + EBV DNA), the proposed nomogram C (including MTBF, nodal necrosis, and the six variables in nomogram B) achieved much better performance in terms of the C-index and ROC analysis. Compared with low-risk patients, high-risk patients stratified by nomogram C had poor DMFS (HR: 4.76, 95% CI: 3.39-6.69; p < 0.0001) and PFS (HR: 4.11, 95% CI: 3.13-5.39; p < 0.0001). With the integration of tumor burden information, nodal necrosis, and other clinical characteristics, MTBF and nodal necrosis were promising supplementary factors to TNM stage and EBV DNA, yielding better performance in predicting survival. As seen in Figure 2A and the coefficients in Supplementary Table 3, MTBF played the most important role in predicting DMFS, and nodal necrosis was also a crucial variable. These results might mainly be associated with two major reasons. First, a high MTBF score and positive nodal necrosis indicate that NPC patients had unfavorable massive primary tumors or regional lymph nodes, which have a greater propensity for occult metastasis (20,31). Second, massive tumors and positive nodal necrosis have been found to be related to the biological aggressiveness of cancer clones, inadequate blood supply, and other adverse factors, including tumor hypoxia (32,33), which is strongly associated with radioresistance and thus relapse and metastasis.
Over 70% of NPC patients are diagnosed with an advanced stage of disease (34,35). According to the National Comprehensive Cancer Network guidelines, locoregionally advanced NPC patients are at high risk for disease progression, for whom ICT + CCRT is recommended as level 2A evidence and CCRT alone as level 2B evidence. However, not all of these patients benefit from ICT (36). In our investigation, our proposed nomogram model C (including MTBF and nodal necrosis) performed well in predicting survival and succeeded in stratifying high-and low-risk patients who benefited from ICT + CCRT and CCRT, respectively. In the high-risk group, the DMFS of patients treated with ICT + CCRT was longer than that of patients treated with CCRT alone (HR 0.60, 95% CI: 0.37-0.97, p = 0.0340; Figure 5A). However, in the low-risk group, the treatment effects of ICT + CCRT and CCRT were similar. This finding might be explained by the fact that ICT not only attenuates tumor load within a brief period to ameliorate tumor hypoxia but also has a systemic cytotoxic effect to eradicate distant micrometastases (37,38). It is well recognized that the improvement in PFS is mainly due to the reduction in distant metastases. In our study, compared with patients who did not receive ICT, patients with high-risk scores treated with ICT showed an improvement in PFS.
Magnetic resonance imaging is now extensively used for the accurate evaluation of NPC to define the TNM stage and follow-up, and the majority of NPC patients undergo MRI before treatment. Our prognostic model utilized the information obtained from pretreatment MRI without increasing the physical or financial burden. One main challenge in our study lies in the fact that three specialists in MRI are required to contour the margin of the primary tumor and regional lymph nodes, which might take time and effort. Fortunately, with the rise of artificial intelligence, deep learning has been widely applied in radiology and pathology for image processing and lesion recognition (39)(40)(41)(42). Automatically recognizing the scope and characteristics of neoplasms will help make our proposed model a more user-friendly prognostic tool. With regard to the result of ROC analysis that the ROC curve of nomogram C vs B was not statistically significant in the two external validation sets, we think the reason lies in the different population distributions and the small sample size. Longer follow-ups and prospective multicenter clinical studies should be carried out to validate our MTBF model and nomogram. Additionally, radiomics has achieved valuable performance in many prediction tasks, including NPC (16,43). In addition to structural tumor burden features, nodal level, and necrosis, further investigation incorporating other radiomics data will be our next research interest. In conclusion, the survival model based on MTBF, nodal necrosis, and clinical factors is a promising prognostic tool for NPC and is helpful for identifying patients who might benefit from ICT.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

ETHICS STATEMENT
This study was approved by the Ethics Committee of the Chinese Clinical Trial Registry (ChiECRCT20190127). Patients' written informed consent was not required because no direct interaction with patients was performed and no personal identification was applied in this study. In addition, this research was conducted in compliance with the Declaration of Helsinki.

AUTHOR CONTRIBUTIONS
XL, CFL, and XiC designed this study. XiC, MQ, GL, JL, MG, ZC, JM, SL, and WX collected data. LK, XuC, KL, YX, and XG rechecked data. XiC, XuC, BJ, and WL performed analysis and wrote the manuscript. CXL, CFL, and XL helped to revise the manuscript. All authors approved the final version of the manuscript.