This study was approved by the Institutional Review Board (IRB) of the First Affiliated Hospital of Zhengzhou University. All patients enrolled in this study provided informed consent.

A total of 115 patients who had a confirmatory diagnosis of esophageal NEC during the period from September 2015 to July 2019 were included in our database. We confirmed the diagnosis of esophageal NEC, without any squamous carcinoma/adenocarcinoma component mixed, on the basis of pathological and immunohistochemical specimen via esophagoscope before surgery or intraoperatively. Patients who had earlier undergone a baseline contrast-enhanced CT before treatment in our institution were enrolled in the study. The remaining 89 patients (median age, 65 years; range, 42-79 years) comprising of 57 males and 32 females were included in the study for treatment observation and subsequent follow-up investigation.

Of the 89 patients, 61 patients (70%) (median age, 65 years; range, 42-78 years) were included in a training cohort for exploring the independent predictors of the texture parameters as well as to optimize the threshold in differentiating the patients’ survival outcomes. The remaining 28 patients (30%) (median age, 65 years; range, 45-79 years) comprised a testing cohort, which was used to validate the prognostic significance of the identified predictors, by applying optimal cutoff values in an external verification set.

The treatment strategies for NEC included esophagectomy (n = 5), chemotherapy (n = 40), radiotherapy (n = 1), chemoradiotherapy (n = 2), preoperative neoadjuvant chemotherapy (n = 1), and postoperative adjuvant chemotherapy (n = 39). The surgical approach for the 5 cases (6%) selected was based on the site of the lesion. The left thoracic approach into the chest was used for a lower esophageal lesion (n = 4), and a right thoracic approach was employed for an upper esophageal lesion (n = 1). The treatment regimen for chemotherapy, preoperative neoadjuvant chemotherapy, and postoperative adjuvant chemotherapy included etoposide (Vepeside, VP - 16) 100mg/m² IV on day 1 to day 5 and cisplatin (DDP) 75mg/m² IV on the first day of a 21 - day cycle for a total of 4 to 6 cycles. Patients treated with radiotherapy received a dose of 1.8 Gray per day on day 1 to day 5, with total dose of 60 Gray at 33 days. Chemoradiotherapy regimen consisted of VP - 16, 50 mg/m² and DDP, 20 mg/m² on days 1 to 5, and a total of 60 Gray of concurrent radiotherapy as mentioned previously.

Therapeutic efficacy was evaluated 12 weeks after surgery by a CT enhancement examination. Patients who were treated with chemotherapy and/or radiotherapy also underwent a CT enhancement examination after every 3 cycles from the beginning of treatment for the first 12 cycles and at every 4 cycles thereafter. All patients underwent clinical follow up through an outpatient service or a telephone communication (every 12 weeks) until death or the end of the study. The OS was set (time from the date of diagnosis to death from any cause) as the main endpoint and was predicted by constructing nomograms.

A baseline CT scan was obtained less than a week before treatment was initiated. All patients fasted for 6 hours. 20 mg of amidoamine was injected intramuscularly 10 - 15 minutes before CT examination. In order to distend the upper gastrointestinal tract, the subjects were first asked to drink 800 - 1000ml water, and subsequently advised to drink 200 - 300 ml water with a straw in a supine position. A contrast-enhanced CT scan of the chest was performed with a 64 - slice CT scanner (Lightspeed VCT, GE Healthcare, Waukesha, WI, USA) or a 256 - slice CT scanner (Revolution CT, GE Healthcare, Waukesha, WI, USA). The scan mode settings comprised of a tube voltage of 120 kV or 80 kV/140 kV with a fast kV - switching technique, tube current under Auto mA with a noise index of 8.0, a slice thickness of 5 mm, a gantry speed of 0.6 seconds per rotation and a pitch of 0.984: 1. A contrast medium containing 350 mg of iodine per ml (Omnipaque™; GE Healthcare, Cork, Ireland) was injected at a flow rate of 3 ml/s via the elbow vein. The dose of the contrast medium was calculated as 1.5 ml per kg body weight. Scanning was triggered when the CT value of the aortic arch reached 100 HU. The venous phase was initiated 60 s after triggering the scan.

The basic morphological image features and clinical data assessed included the age and gender of the patient, tumor length, location, the thickness of esophageal tumor wall (Median), clinical symptoms, immunohistochemical index of markers such as Cytokeratin (CK), Chromogranin A (CgA), Synaptophysin (Syn), Neuron-specific enolase (NSE), and Ki-67 (%). Tumor invasion depth (cT stage), lymph node metastasis (cN stage), distant metastasis (cM stage), treatment strategy, tumor margin (well-defined or ill-defined), presence or absence of tumor calcification, morphological subtype and tumor homogeneity (homogeneous or presence of necrosis) were also assessed. The degree of enhancement was classified as slight, moderate, and marked enhancement, with the standard used for assessment of enhancement being: slight enhanced, the enhancement of the nodule was close to that of adjacent muscles; moderate enhanced, enhancement slightly higher than the adjacent muscles; marked enhanced, enhancement markedly higher than in the muscles (14). The start date of treatment, and the date of pretreatment CT examination were also documented.

A radiologist (Y.Z.) with 12 years of experience in chest and abdominal CT diagnosis, independently retrieved the CT images in a DICOM format from the picture archiving and communication system (PACS), and resampled them to 1 × 1 × 1 mm³ sizes for maintaining a consistent thickness. The resampling images were subsequently loaded into the ITK-SNAP software (version 3.6.0; www.itksnap.org) for manual segmentation. Tumor volume of interest (VOI) segmentation was performed, and the images were exported in an NII. format for subsequent 3D reconstruction, data extraction, and analysis. First, the VOI was delineated to encompass the whole enhanced esophageal lesion in the sagittal position. Visually identified air density within the tumor and the normal esophageal wall adjacent to the margin of the lesion were excluded. Next, the outlined segmentation was compared and modified in an axial position and coronal position, respectively. A 3D reconstruction was precisely performed along the edge of the VOI region which was automatically separated from the other structures and the background.

The textural features were derived from a multiparadigm numerical computing and programming software Artificial Intelligence Kit (version 3.2.2, GE Healthcare). A total of 402 radiomics features were extracted for tumor characteristics from the VOI region segmented from the pretreatment venous phase CT images. The features were composed of the histogram (42 features), gray level co-occurrence matrix (GLCM) (144 features), gray level size zone matrix (GLSZM) (11 features), gray level run length matrix (GLRLM) (180 features), formfactor (15 features), and Haralick (10 features) (16–21). An analysis of variance (ANOVA) was performed for dimension reduction, and a correlation analysis was performed to reduce data redundancy. The details of features extraction are described in the supplementary data . The extracted texture features were standardized so as to remove the unit limits of the data of each feature.

A random selection of 30 enhanced CT images of the patients was selected for assessment of intraobserver agreement (20). The intraclass correlation coefficients (ICCs) were calculated to evaluate inter-observer variability of radiomics features and were classified as poor (< 0.40), fair (0.40 - 0.59), good (0.60 - 0.74), or excellent (0.75-1.00) (22). Highly correlated features (ICC ≥ 0.90) were identified and selected initially. After a spearman correlation analysis between the texture features, there were a remainder of 50 features.

Potential prognostic factors of survival were evaluated using a univariate Cox proportional hazard model based on forward stepwise selection from the clinical variables. A least absolute shrinkage and selection operator (LASSO) Cox regression method was used to determine features that predicted the overall survival. The radiomics score for each patient was computed using an equation in which the coefficients were derived from the LASSO Cox model. The radiomics score (Rad-score) was computed for each patient through a linear combination of selected features weighted by their respective coefficients. Feature selection and the radiomics signature building were performed in training and testing sets. A five-fold cross-validation test was used to detect the goodness of fit of the model. The model was validated for its calibration ability by calculating the probability of outcome for each patient of the whole dataset according to the model and comparing it with the actual survival of the patient. More details of LASSO Cox model and radiomics signature building are described in the Supplementary Data . The survival models used different sets of factors that included the radiomics score only, and a combined score using both selected radiomics features and clinical variables.

Statistical analysis was conducted with R software (version 3.4.4; http://www.Rproject.org) and p < 0.05 was considered statistically significant. The differences between the continuous variables were compared using the Mann-Whitney U test, and the differences between the categorical variables were analyzed by the Chi-square test. The predictive accuracy of the radiomics signature was quantified by receiver operator characteristic (ROC) curve in both, the training and the testing sets. The median value of each score in the model was used to stratify the patients into high-risk and low-risk groups (23). The potential correlation of the two models with the overall survival (OS) was assessed using Kaplan-Meier analysis. The difference in survival between the groups was calculated as hazard ratio (HR) for each model by use of log-rank tests.

The details of patient enrollment are summarized in ( Figure 1 ). The basic clinical characteristics of patients in the training and testing cohorts are shown in ( Table 1 ). There was no statistical difference between the two cohorts in terms of the variables gender, age, morphological subtype, tumor length, location, clinical symptom, immunohistochemical result, the thickness of esophageal tumor wall, Ki-67 (%), cTNM stage, treatment strategy, tumor margin, necrosis, calcification, enhanced homogeneity, enhanced pattern, and enhanced degree. The median OS in the training and testing cohorts were 41.7 months (range: 33.4 - 50 months) and 45 months (range: 24.6 - 65.4 months), respectively (p = 0.302). OS was censored in 18 and 6 patients from the training and testing cohorts, respectively.

The typical example of the drawn VOI covering the largest image slice of the primary lesion is depicted in ( Figure 2 ). The inter-observer reproducibility analysis was measured by ICC on the VOI-based feature extraction, with a mean ICC of 0.67 (range: -0.03 to 1). The LASSO regression was used for data dimension reduction, feature selection, and radiomics signature building. After feature selection derived from the LASSO regression, the COX regression clinical prediction model was constructed, and the nomogram was created to predict individuals. The result of data dimension reduction with LASSO Cox regression models is described in ( Figures 3A, B ). The seven selected radiomics features consisted of mean deviation, GLCM Entropy_angle90_offset1, inverse difference moment, high Grey Level Run Emphasis_All Direction_offset4_SD, Run Length Nonuniformity_All Direction_offset4_SD, Short Run Emphasis_All Direction_offset1, and Sphericity. The values of the seven selected features in each patient were applied into the formula, and the Rad-scores were obtained to reflect the risk. The Rad-scores thus derived by a combination of the seven most valuable variables and their coefficients are shown in ( Figure 3C ).

For clinical variables, as described in ( Table 2 ), the univariate Cox regression analysis indicated that only the enhancement degree was an independent factor that influenced the OS. For the combined model with both radiomics features and clinical variables, the combined score was calculated using selected features and variables. The nomograms for radiomics-based model and the combined model were built based on the training set. One-year, 3-year and 5-year survival probabilities estimated with the nomograms, indicated a good fit for both models in the training and testing sets as seen in ( Figure 4 ).

Patients were divided into high and low risk groups using the predicted risk of the model. For the radiomics based model and the combined model, the thresholds were adopted as the median values of the distribution of the risk groups for OS in both, training and testing datasets. To find the best threshold to make the best classification in radiomics based model, the highest value of the Youden index is selected. Threshold values were -0.234 for the radiomics based model. As depicted in ( Table 3 ), the radiomics signature based model showed a favorable predictive efficacy than the combined model that integrated the radiomics signature with clinical risk factor, with a C-index of 0.844 (95% CI: 0.783 - 0.905) and 0.847 (95% CI: 0.782 - 0.912) in the training set, and a C-index of 0.805 (95% CI: 0.707 - 0.903) and 0.745 (95% CI: 0.639 - 0.851) in the testing set, respectively. The difference was statistically significant in the training and testing datasets (p < 0.0001 for each comparison). The corresponding Kaplan-Meier curves are shown in ( Figure 5 ).

This study has several limitations. Firstly, this is a retrospective study with a small sample size, which might have caused instability in features selection in both, the training and the testing sets. Future studies with larger samples are needed to thoroughly validate the results of this study. Secondly, the results were derived from the use of two CT scanners with a different number of detectors, which might have contributed to a bias in radiomics features extraction. Moreover, additional imaging modalities, such as radiomics features derived from pretreatment functional MR images can provide information on different aspects of tumor characteristics that can be crucial for predicting tumor prognosis. Future research in this context will aim at the correlation between the imaging features and gene or protein biomarkers for survival prediction.