All protocols of this prospective and observational study were approved by the Institutional Review Board (Project ID: 201404005). Patients from Shandong Cancer Hospital and Institute were randomly assigned to the training set and the internal testing set in a ratio of 7:3 to ensure an adequate sample size. The external testing set was enrolled in Anyang Tumor Hospital.

All methods were carried out following relevant guidelines and regulations. All patients provided written informed consent before enrollment. Inclusion criteria were as follows: (a) ESCC histologically proven by endoscopic biopsy; (b) no contraindications to MR examination; (c) clinical stage T3 or T4 defined on endoscopic ultrasonography, diagnostic CT scans, or ¹⁸fluorodeoxyglucose positron emission tomography/computed tomography (¹⁸F-FDG PET/CT) scan; (d) no prior anticancer treatment; (e) no distant metastasis except lymph nodes; (f) no other primary tumor. Clinical staging was based on the tumor-node-metastasis classification of malignant tumors (UICC, Version 8th). All patients were treated with concurrent CRT. The daily fractional dose of radiotherapy was 1.8–2.0 Gy, administered 5 days a week, and the total dose was 50.4–66.0 Gy with 6–15 MV X-rays performed by intensity-modulated radiation therapy. Before the early treatment response evaluated, all patients accepted two cycles of chemotherapy concurrently during radiotherapy every 3 weeks. One of the following proposals of chemotherapy would be adopted: 1) an intravenous injection of paclitaxel (150 mg/m²) on day 1 and cisplatin (75 mg/m²) on days 1–3; and 2) S-1 administered orally twice daily at 80 mg/m² for 2 weeks and cisplatin (75 mg/m²) on days 1–3.

The response was assessed according to RECIST 1.1 (6). Upper gastrointestinal endoscopy, chest and abdomen CT, or whole-body ¹⁸F-FDG PET/CT was performed 2–3 months after treatment finished to compare with the pretreatment images. The definition of chemoradiotherapy sensitivity was as follows: the sensitive group comprised patients with complete response (CR) and partial response (PR); the resistant group consisted of patients with stable disease (SD) and progression disease (PD).

MRI examination with DWI and T2-weighted imaging (T2WI) in the axial orientation was sequentially performed at about 1–3 days before treatment, 5^th fraction of radiation completed (5f), and 10^th fraction of radiation completed (10f) as shown in Figure 1A . All patients were scanned using a 3.0 T MRI system (Achieva, Philips Medical Systems, Netherlands; Ingenia, Philips Medical Systems, Netherlands) with body phased-array coil anterior and spine array coil posterior. T2WI was performed by using a two-dimensional fluid-attenuated inversion recovery sequence (repetition time/echo time, 1,277 ms/70 ms; flip angle, 140°; bandwidth, 315 kHz; section thickness, 4 mm; slice gap, 1 mm; field of view, 350 × 350mm; matrix, 384 × 276; no. of slices, 24–36; imaging time, 5 minutes). DWI was performed by using a single-shot echo-planar imaging (repetition time/echo time, 1,244.7 ms/47.9 ms; flip angle, 140°; bandwidth, 260 kHz; section thickness, 4 mm; slice gap, 1 mm; field of view, 360 × 360 to 400 × 400 mm; matrix, 128 × 128; no. of slices, 24–36; imaging time, 5 minutes). The number of signal averaged (NSA) was 1. Diffusion gradients were applied in three orthogonal directions by different diffusion weightings (b-values) of 0, 300, and 600 s/mm². Especially, b = 600 s/mm² was assigned a NSA of 4 to improve the signal-to-noise ratio. Respiratory triggering was used.

The ADC map was generated automatically at the MR workspace (Philips Medical Systems Extended MR workspace, Netherlands) from DWI, which was overlayed and averaged in three directions with the b-value of 0 and 600 s/mm².

The location of the primary tumor was firstly identified on DWI (b = 600 s/mm²) as areas of high signal ( Figure 2 ). Although it was not outlined on DWI, it assisted in accurately finding the corresponding location of the tumor in the ADC map. Regions of interest (ROIs) were manually drawn on the ADC map via a segmentation editor in the 3D Slicer software (version 4.10.0, http://www.slicer.org) (23). The definition of ROI was the region of the tumor that was of low signal in the ADC map but excluding the lumen. Every slice containing the tumor was outlined. These ROIs constituted a volume of interest (VOI) with the entire tumor visualized on the ADC map as showed in Figure 2 . The mean tumor ADC value was calculated by averaging the measured ADC values of the VOI. For whole-tumor radiomics feature, the Slicer Radiomics, which was based on PyRadiomics program (Revision 2.2.0) in the 3D Slicer software, automatically extracted whole-tumor RFs from the three-dimensional reconstructed VOI (24). In our practice, resampled voxel size, LoG kernel sizes, and bin width were respectively set to 1, 1, and 25. In total, 851 RFs were extracted including 18 first-order statistics features, 14 shape-based (2D and 3D) features, 24 gray-level co-occurrence matrix (GLCM) features, 16 gray-level run length matrix (GLRLM) features, 16 gray-level size zone matrix (GLSZM) features, 5 neighboring gray tone difference matrix (NGTDM) features, 14 gray-level dependence matrix (GLDM) features, and 744 wavelet-based features. Detailed classification of the extracted RFs was summarized in Supplementary Table 1 . In addition, the RFs mentioned in our study were consistent with the definition of the image biomarker standardization initiative (25).

Selection procedures were performed on the pretreatment images with three steps.

Because DWI is susceptible to field inhomogeneities, large amounts of magnetization artifacts in the chest area, and patient movement, two consecutive DWI acquisitions of the same patient may give images with slightly different imaging characteristics, which may affect radiomics analysis. Thus, the first is to evaluate the repeatability of the radiomics features between two ADC maps in a short time by the intraclass correlation coefficient (ICC) method. Five extra patients with ESCC were enrolled before the first treatment in Shandong Cancer Hospital and Institute ( Supplementary Table 2 ) who were not included in the training and internal testing sets. The interval between the two consecutive DWI scans was 5 minutes ( Supplementary Figure 1 ). VOIs were delineated in the ADC map by an experienced radiologist (WL, with 20 years of experience in MRI), then the VOI was duplicated to the other ADC map from the same patient. The reliability coefficient of RFs was supposed to be higher than 0.75.

Secondly, to select RFs with high interobserver reproducibility, VOIs were delineated by two experienced radiologists (WL, appointed as reader 1; NS, appointed as reader 2 with 11 years of experience in MRI) who were blind to clinical information, treatment plan, and response about patients. Bland–Altman analysis was used to test the interobserver reproducibility of the RFs from the VOI delineated by different radiologists. The differences between the same parameters extracted from two readers’ delineation were plotted against their average and reported as a percentage. The lower and upper reproducibility limits were calculated as ± 1.96 standard deviations. In our case, the mean and standard deviation of the differences were supposed to be less than 5% and 10%, respectively, as stated in related work (26, 27).

Thirdly, to reduce the number of features and select the most significant features correlated with treatment response, the minimum redundancy maximum relevance (mRMR) algorithm was used to identify and rank the top 30 features. Finally, the selected RFs were applied in the training, internal, and external testing sets from the delineation of reader 1.

The changes of the whole-tumor ADC value between images of different time points were respectively calculated as follows:

where ΔADC_{1st week}, ΔADC_{2nd week}, and ΔADC_{2 weeks} denoted the change of the ADC value within the 1st week, the 2nd week, and the first two weeks, respectively. Similarly, ADC_pretreatment, ADC_5f, and ADC_10f denoted the ADC value of images acquired before treatment, after the 5th radiation, and after the 10th radiation, respectively.

The changes of the selected RFs were defined as the relative difference between each measurement as follows:

where ΔRF_{1st week}, ΔRF_{2nd week}, and ΔRF_{2 weeks} denoted the change of RFs within the 1st week, 2nd week, and first two weeks, and RF_pretreatment, RF_5f, and RF_10f denoted the RFs in images acquired before treatment, after the 5th radiation, and after the 10th radiation, respectively.

As shown in Figure 1B , RF_pretreatment, ΔRF_{1st week}, ΔRF_{2nd week}, and ΔRF_{2 weeks} were used to train a support vector machine (SVM) classifier in the training set, and then a radiomics signature was built by the trained classifier and evaluated in the internal and external testing sets. The cost‐based SVM classifier, which used the radial basis function kernel and a cost coefficient of 0.001, was trained by RFs with the 10-fold cross-validation. The SVM classifier transformed the feature space to a high-dimensional space where a separating hyperplane maximized the margin between classes. For each patient, the SVM classifier generated a radiomics signature for evaluating predictions of treatment response. The receiver operating characteristic (ROC) curves of obtained radiomics signatures in the training, internal, and external testing sets were analyzed and then tested by the DeLong test. The chi-square test was used to compare the predicted responses, which were corresponding to radiomics signatures in the internal testing with the signatures’ ground truth, i.e., treatment response evaluated by RECIST1.1. Afterward, the eligible radiomics signatures together with clinical characteristics were imported to a multivariate logistic regression-based prediction model to predict the sensitivity to treatment. Furthermore, the performance of radiomics signatures, the ADC value, clinical characteristics, and the regression model was assessed by the ROC curve, respectively. The area under the curve (AUC) was calculated, and the optimal Youden’s index was determined.

Statistical analysis was performed using the R software (version 3.6.1, http://www.R-project.org). The ICC method, Bland–Altman method, mRMR approach, SVM classifier, logistic regression analysis, and ROC curve analysis were implemented based on the R packages “irr”, “BlandAltmanLeh”, “mRMRe”, “e1071”, and “pROC.” Chi-square test or t-test was used to examine the correlations between clinical features and treatment response. The differences in ADC and ΔADC values were compared between the sensitive group and the resistant group by Mann–Whitney U test. P-value < 0.05 was considered as an indicator of the statistically significant difference.

In Shandong Cancer Hospital and Institute, eighty-five consecutive patients with ESCC were enrolled. Nine patients were excluded, given that the DWI of 6 patients were deformed or had high noise and 3 patients’ radiotherapy was interrupted. Finally, 76 patients were enrolled in the study during the period of enrollment from June 2014 to September 2019. Fifty-four of the 76 patients experienced clinical PR, 2 with CR and 20 with SD. Patients (n=76) were randomly assigned to training (n=53) and internal testing sets (n=23) in a ratio of 7 to 3. There was no significant difference in clinical characteristics between the training and internal testing sets (sex, age, tumor location, clinical T stage, lymph node status, and radiation dose) and treatment response in both sets ( Table 1 ). In Anyang Tumor Hospital, twenty consecutive patients with ESCC were enrolled during the period of enrollment from April 2017 to May 2019. Three patients were excluded, given that 3 patients’ radiotherapy was interrupted. Finally, 17 patients were enrolled in the external testing set. Eleven of the 17 patients experienced clinical PR and 6 with SD. Clinical characteristics are presented in Table 1 .

We measured the voxel volume of the whole-tumor ADC map to substitute the tumor volume and the maximum 3D diameter and their changes between each time point. Supplementary Tables 3 , 4 show that these parameters had no relationship with treatment response.

The mean values of ADC_pretreatment, ADC_5f, and ADC_10f of all patients are shown in Table 2 . All of the ADC values and relative changes during the first two weeks showed no significant difference between the sensitive and resistant groups in the training and internal testing sets. In the external testing set, ADC in the 5^th radiation showed association with treatment response (p = 0.048). We noticed that the ADC value was gradually higher than before treatment as the treatment progressed.

A total of 560 RFs ( Supplementary Table 5 ) were selected by the ICC method from two consecutive whole-tumor ADC map acquisitions in the same patient. Afterward, 224 RFs ( Supplementary Table 6 ) were selected from pretreatment images in the training set via the Bland–Altman method, and 30 RFs ( Supplementary Table 7 ) were finally chosen out according to mRMR. Most of these RFs were calculated to measure the local homogeneity of the image. To predict good and bad responder, four radiomic signatures based on an SVM classifier were respectively built with RF_pretreatment, ΔRF_{1st week}, ΔRF_{2nd week}, and ΔRF_{2 weeks}. The effectiveness of these 4 radiomics signatures to classify the sensitive group versus the resistant group is shown in Figure 3 . Only the radiomics signature based on ΔRF_{2 weeks}, denoted as R-Signature_{2 weeks}, discriminated treatment response with an AUC value higher than 0.5 in the training set [AUC = 0.824, 95% confidence interval (CI): 0.679–0.968], the internal testing set (AUC = 0.744, 95% CI: 0.465–1.0), and the external testing set (AUC = 0.742, 95% CI: 0.478–0.919). No difference was found between the training and internal testing sets according to the results of the DeLong test (p = 0.580), nor between the training and external testing sets (p = 0.526). Compared with the evaluated response, the performance of response prediction in internal testing was acceptable (p = 0.048). The chi-square test of prediction results between R-Signature_{2 weeks} in the training set and the internal testing set achieved a p-value of 1.000, which implied that the distribution of R-Signature_{2 weeks} between the two sets was not statistically significantly different.

As elaborated in Supplementary Table 8 , in the training set, the outputs of univariate logistics regression showed that tumor location and R-Signature_{2 weeks} were separately associated with treatment response. We performed the tumor location and R-signature_{2 weeks} to the prediction model of treatment response by multivariate logistics regression ( Supplementary Table 9 ) in the three sets. R-Signature_{2 weeks} was identified as an independent factor predicting treatment response in multivariate analysis in the three sets. Figure 4 shows the results of ROC analyses on the tumor location, R-Signature_{2 weeks}, and the regression model in the prediction of treatment response in the training, internal, and external testing sets. Table 3 shows the effectiveness of our proposed regression model in discriminating the sensitive and resistant groups.