Computed tomography-based radiomic model predicts radiological response following stereotactic body radiation therapy in early-stage non-small-cell lung cancer and pulmonary oligo-metastases

1. Patient selection

Patients receiving thoracic SBRT using active breathing control (ABC) technique between January 2008 to December 2018 were retrospectively analyzed. The study obtained the approval of the Institutional Review Board of the University of Hong Kong/Hospital Authority Hong Kong West Cluster (HKU/HA HKW) (No. UW 16-522).

Inclusion criteria were stage I NSCLC (i.e., the American Joint Committee on Cancer staged T1-2N0M0) and pulmonary oligo-metastases which had primary tumors surgically removed, had 5 or fewer lung metastases and stable systemic control.

Exclusion criteria were inability to tolerate the breath-hold required for ABC technique, and thoracic tumors within the no-fly zone as defined by Timmerman et al. [11] or tumors which overlapped with mediastinal/chest wall structures. The above criteria were chosen to avoid confounding factors which make interpretation of radiomic textual features inaccurate, such as tumors with excessive movement, tumors treated with alternative SBRT dose schedules and tumors overlapping with other structures.

A total of 69 patients were recruited. All patients received 45–60 Gy in 3–5 fractions on alternate days equivalent to BED₁₀ >100 Gy. Clinical parameters of the patients were summarized in Table 1.

Table 1.

Patient and treatment characteristics

Patients’ breathings were actively controlled using ABC that employed a modified spirometer and occlusion valves. By closing both valves at predefined lung volume, breathing was held at a definite and stable expansion for 20 seconds while radiation beam was on. Patient was then allowed to breathe freely until the next ABC cycle. CT simulation was performed using 3 mm slice thickness.

Treatment planning was carried out with the Eclipse radiotherapy planning system (Varian Medical Systems, Palo Alto, CA, USA). The gross tumor volume (GTV) was delineated using lung windows at width of 1,600 and level of -600. The planning target volume (PTV) was defined as GTV plus 5 mm margin. A 10–20 Gy per fraction was prescribed at 80%–90% isodose level of PTV to a total of 45–60 Gy in 3–5 fractions over 2 weeks. Prescribed dose had to cover 95% of PTV and 99% of PTV received >90% of prescribed dose. Treatment isocentre was verified before each fraction and compared to center of the tumor on CT image.

2. Outcome assessment

The primary outcome of the study was local radiological response, classified as radiological CR, partial response (PR), stable disease (SD), and progressive disease (PD).

Other outcomes such as locoregional failure, distant metastasis and overall survival were also studied. Survival time was calculated from time of initial irradiation.

Regular thorax CT with contrast at 3 months interval in the first year and at 6 months interval from second year onward were performed for all patients in our tertiary institute. Additionally, PET/CT could be performed as per clinician’s discretion. Tumor response was assessed by independent senior radiologists according to the Response Evaluation Criteria in Solid Tumors (RECIST) version 1.1 and Green’s criteria [12]. Tumors with stable or progressive disease after SBRT were grouped together as non-responders while tumors with CR and PR were grouped together as responders. Each category of response was analyzed using the best response achieved throughout the entire follow-up period.

Progression of disease within the PTV was defined as local recurrence. Locoregional recurrence was defined as recurrence in the lungs and or in regional lymph nodes. Distant metastasis was specifically defined according to the stage at diagnosis. For early-stage lung cancers, it was defined as development of lesions to distant organs or lymph nodes. For metastatic malignancies, it was defined as development of new metastases in addition to the original extent of disease before SBRT. Patients were censored regularly at least every 3 months and also when new follow-up imaging was available.

3. Radiomic and statistical analysis

The contoured CTs were analyzed using 3D-Slicer software. The voxel spacing of the CTs was normalized to 1 mm × 1 mm × 1 mm. Bin width of 25 HU (Hounsfield unit) was applied to the images for radiomics feature extraction.

One hundred and ten radiomic features of the GTVs including shape, statistical and textual features were extracted using PyRadiomics module [13] (Supplementary Table S1). The extracted features were subsequently analyzed using in-house software developed in MATLAB R2020a (MathWorks Inc., Natick, MA, USA). To determine features predictive of non-responder, the receiver operating characteristic (ROC) curve was computed for individual radiomics features to distinguish non-responders (SD and PD) and responders (PR and CR). The area under curve (AUC) was calculated for each ROC curve. The significance of AUC values were evaluated using z test with adjustment for multiple testing by Bonferroni adjustment [14]. Cox regression were also done to elucidate the clinical factors that affects overall survival.

A model using support vector machine (SVM) was trained to predict treatment response based solely on the extracted radiomics features in accordance with the TRIPOD statement [15]. The model aimed to classify patients’ radiological response into three classes, namely non-responder (PD and SD), PR, and CR.

In the present study, patients with NR outweighed those with PR and CR. As the response grouping was unbalanced, adaptive synthetic sampling method (ADASYN) was used to up-sample minority classes (PR and CR) to improve model performance. All extracted radiomics features were included in the training of the SVM model without prior selection to prevent bias. No clinical parameters were included in order to assess the ability of the model to predict response solely based on radiomic features.

Gaussian SVM was used in this study with one versus one classification technique and lasso regularization. Hyperparameters such as box constraint levels and kernel scale were optimized using Bayesian optimization during training. The hyperparameters of the best performing model would be selected. As radiomic features are highly correlated and at time redundant, lasso regularization was employed. During the process of training, weighting coefficients of different features were computed with lasso regularization. As a result, highly important features would have greater contribution towards the final model while redundant and noisy features would have little if any influence on the final model. Overfitting of SVM model can limit the generalizability of model to external datasets. Hence, internal validation through 10-fold cross-validation was used to avoid overfitting while maintaining accuracy. ROC curves were constructed with bootstrap resampling for 1,000 times. AUC was used to evaluate model performance. Schematic representation of model training was presented in Fig. 1.

Fig. 1.

Radiomic workflow in the present study. Radiomic features (including shape, statistical, and textual features) were extracted from contoured tumor using PyRadiomics module. Under-represented population (CR and PR) was up-sampled using ADASYN algorithm. A Gaussian SVM model was then trained with lasso regularization. Hyperparameters were tuned using Bayesian optimization. Internal validation was done by 10-fold cross-validation. The performance of the model was then evaluated and prediction generated. CR, complete response; PR, partial response; ADASYN, adaptive synthetic sampling method; SVM, support vector machine.

Statistical analyses were conducted by Statistical Package for the Social Sciences (SPSS) version 25.0 (IBM, Armonk, NY, USA). Statistical significance was tested using t-test for continuous variable and chi-squared test for categorical variables. Log rank test was used for univariate analysis. A p-value less than 0.05 was considered significant.

1. Clinical characteristics

Sixty-nine patients were recruited in total from 2008 to 2018. Median age was 73 years (range, 44 to 92 years). Forty-two patients were male while 27 patients were female. Out of the 69 patients included, 55 patients had lung primary while 14 patients had oligo-metastases to the thorax from other tumors. For the 55 patients with lung primary, 29 patients had non-metastatic lung cancer, while 26 patients had metastatic lung cancer. Overall, 40 patients had metastatic disease while 29 patients had non-metastatic disease.

As some patients had more than one irradiated tumor, a total of 85 lesions were analyzed, of which 31 achieved CR, 11 achieved PR, and 43 tumors were non-responders (NR). As 43 tumors out of 85 lesions were non-responders, the objective response rate was lower than that reported in previous studies [11]. Most extra-thoracic malignancies was colorectal in origin (16 tumors). Other extra-thoracic diseases included hepatocellular carcinoma, head and neck and esophageal cancers. In terms of histology, the large majority was represented by adenocarcinoma (72 tumors), while others were squamous cell carcinomas (n = 7), hepatocellular carcinomas (n = 5), and small cell carcinoma (n = 1).

Median follow-up time was 50.5 months (range, 1 to 132 months). Median overall survival post-SBRT was 32 months (range, 1 to 132 months; confidence interval [CI], 24.72–39.28).

Median local failure-free survival (LFS) was 31 months (range, 1 to 132 months; CI, 23.6–38.4). Median locoregional failure-free survival (LRFS) was 22 months (range, 1 to 132 months; CI, 14.46–29.54). Three-year local control rate was 87.15%, while 5-year local control rate was 80.12%. Median distant metastasis-free survival (DMFS) was 18 months (range, 1 to 132 months; CI, 10.46–25.54).

For NSCLC patient, overall survival was 38 months (range, 1 to 132 months; CI, 29.7–46.3), LFS was 37 months (range, 1 to 132 months; CI, 29.3–45.9), LRFS was 36 months (range, 1 to 132 months; CI, 28.1–45.2), and DMFS was 35 months (range, 1 to 132 months; CI, 27.1–44.2). Further subdividing patient into metastatic lung cancer and non-metastatic lung cancer. For metastatic lung cancer patients, OS was 45 months (range, 9 to 116 months; CI, 33–65), LFS was 36 months (range, 9 to 116 months; CI, 29–62), LRFS was 36 months (range, 9 to 116 months; CI, 29–62), and DMFS was 36 months (range, 9 to 116 months; CI, 29–62). For non-metastatic lung cancer patients, OS was 33 months (range, 1 to 132 months; CI, 35–62), LFS was 33 months (range, 1 to 132 months; CI, 35–62), LRFS was 33 months (1 to 132 months; CI, 35–62), and DMFS was 33 months (range, 1 to 132 months; CI, 35–62).

Regarding non-NSCLC patients, OS was 48 months (range, 13 to 116 months; CI, 32.6–67.3), LFS was 36 months (range, 9 to 116 months; CI, 27.8–63.7), LRFS was 36 months (range, 9 to 116 months; CI, 27.8–63.7), and DMFS was 30 months (range, 13 to 116 months; CI, 14.4–45.6).

The median overall survival for responders and non-responders to SBRT was 47 months (range, 9 to 132 months; CI, 36.3–57.3) and 33 months (range, 1 to 93 months; CI, 23.93–42.07), respectively (Fig. 2) with trend towards improved survival for responders (CR and PR) versus non-responders to SBRT (p = 0.059).

Fig. 2.

Kaplan-Meier curve demonstrating trend towards improved overall survival in responders compared to non-responders to stereotactic body radiation therapy.

The trend towards improved overall survival in responders versus non-responders to SBRT (p = 0.062) was confirmed by multivariate Cox regression analysis. No significant relationship with overall survival were demonstrated for other clinical characteristics (Table 2). This may be due to the small sample size in the present study. All treated tumors received SBRT dose ranging from 45–60 Gy in 3–5 fractions, equivalent to BED₁₀ >100 Gy. SBRT dose across non-responding and responding groups were similar. Patient and tumor characteristics were summarized in Table 1. All clinical characteristics including patient, disease and treatment features were not significantly different across groups.

Table 2.

Cox regression for overall survival

2. Radiomics features

We first sought to determine radiomic features that can predict response to SBRT by univariate analysis. Tumors were divided into non-responders (SD and PD) and responders (CR and PR). ROC curves were constructed to predict tumor response.

Out of the 110 radiomic features, skewness and root-mean-square (RMS) were identified as statistically significant markers of NR to SBRT (AUC, 0.619–0.629; p < 0.05) (Table 3). Both features are first order statistical features describing the statistical distribution of voxel intensities. Skewness is a statistical features describing the skewed distribution of voxel intensities with the region of interest. Positive skewness indicates an elongated tail to the right side of the mean on voxel intensity histogram. RMS refers to the RMS of voxel intensities within the region-of-interest (Table 3). All clinical parameters were non-predictive of treatment response (p > 0.05).

Table 3.

Significant radiomic features in univariate analysis

3. Multivariate machine learning model

As the predictive accuracy of single radiomic parameter was unsatisfactory, contemporary multivariate machine learning models may be able to improve prediction accuracy. After up-sampling the minority classes (PR and CR) using ADASYN, the SVM model was trained to classify NR, PR, and CR using one versus all technique. The overall accuracy of classification was 74.81%. F1 scores for CR, PR, and NR were 0.71, 0.84, and 0.72, respectively. The AUC for CR, PR, and NR prediction was 0.86 (95% CI, 0.794–0.921; p < 0.05), 0.946 (95% CI, 0.873–0.978; p < 0.05), and 0.857 (95% CI, 0.789–0.915; p < 0.05), respectively (Fig. 3).

Fig. 3.

Receiver operating characteristic curve for support vector machine prediction of radiological response: (A) complete response, (B) partial response, (C) non-responders. AUC, area under curve.

1. Use of radiomics to improve prognostication in SBRT

With the advent of precision medicine, more accurate clinical models have been devised for prognostication in oncology. Most tumors are highly heterogeneous in nature in terms of its genomic expression, associated stroma and vasculature and this has profound implication on prognosis [16]. Variability in not only inter-individual, but as well remarkable among different lesions within the same individual and in different areas of the same lesion, as in metastatic tumors [17]. Hence, clinical features may not completely capture the heterogeneity of tumors.

The advent of radiomics might partially overcome this limit, as clinically undetectable features can be extracted to capture intra-tumoral and across tumor heterogeneity [10]. Hence, efforts have been made to harness radiomics for prognostication.

Huynh et al. [18] demonstrated that radiomic features can be prognostic for overall survival and distant metastasis in patients with early-stage lung cancer receiving SBRT. Multi-objective models such as random forest were constructed subsequently to improve prognostication in early stage lung cancer [19].

Apart from lung cancer, radiomic driven models have been applied in various other malignancies such as colon, prostate and breast cancer [20-22]. Both treatment outcomes and complications have been accurately predicted for various primary sites [23-25]. Specific to SBRT, clinical outcomes such as overall survival, distant metastasis or complications such as pneumonitis have been modelled using radiomics [26,27]. It was also applied in other imaging modalities such as PET/CT or MRI [28].

SBRT is commonly employed for early-stage lung cancer with curative intent, and is increasingly used in patients with metastatic lung cancer with oligo-progressive or oligo-metastatic disease, allowing high local control rates [29]. This approach can achieve superior local control compared to conventional radiotherapy [30,31]. However, a pathological series by Palma et al. [4] revealed a 60% pCR rate after SBRT for early-stage NSCLC, which is lower than previously hypothesized. The rate of local control may be lower in other unfavorable histology or in more advanced disease [32]. The choice of surgical treatment versus SBRT is at the discretion of the clinician and should be based on patient’s characteristics, comorbidities and on clinical experience, preferentially in the context of a multi-disciplinary team [33].

Therefore, pre-SBRT treatment outcome prediction has pragmatic implication in the optimal selection of treatment modalities. Recent radiomic models allowed accurate prediction of overall survival, DMFS, and progression-free survival for early-stage NSCLC receiving SBRT [18,34].

The population of this study included patients with both early and metastatic disease as well as different primary site and histology, with the aim is to identify universal radiomic markers of radiological response to thoracic SBRT regardless of primary site and histology.

It is of note that, despite including both metastatic and non-metastatic malignancies, there was a trend towards improved survival for radiological responders compared to non-responders. As such, radiological response may be an important prognostic marker for thoracic SBRT that could reflect an impact on overall survival. Thus, anticipating radiological response prior to administrating SBRT could play a relevant role in treatment decision.

2. Predictive radiomics features

To our knowledge, this is the first study recruiting patients receiving SBRT including both primary lung cancer and lung metastases from extra-thoracic disease. Skewness and RMS were identified as predictors of radiological non-responders. Interestingly, these two features were also identified in previous radiomic studies on conventional thoracic radiotherapy. Coroller et al. [35] identified seven features that were significantly associated with pathological gross residual disease following chemo-irradiation in early-stage lung cancer including skewness and RMS. These findings were confirmed by a study by Chong et al. [36] demonstrating that skewness was predictive of pathological non-responders following chemo-irradiation in lung adenocarcinoma. Considering other primary tumors, skewness was as well found to be associated with poorer prognosis in patients with colorectal cancer [37].

Changes in skewness over time were also found to be informative; longitudinal decrease in skewness post treatment was demonstrated to be predictive of better overall survival in lung adenocarcinoma [38]. As increased skewness has been found to confer poor prognosis in multiple radiomic studies across multiple tumor sites and treatment modalities, it is possible that increased skewness could be a universal negative radiomic prognostic factor. The possible underlying biological mechanisms are further discussed in the subsequent section.

Although less data exists regarding the prognostic implication of RMS, it was recently demonstrated that lower RMS was associated with poorer survival in lung cancer patients [39].

Despite enrolling a small cohort, the present study is in accordance with previous radiomic studies. In addition, the findings were confirmed in different studies utilizing different imaging and analysis protocols [35,36,39]. This suggests that skewness and RMS are potentially important poor radiomic prognostic factor with generalizability.

3. Biological basis of radiomic features: radiogenomics

Increasing efforts have been made to uncover the biological basis underlying radiomic features, as biological alterations may have specific radiological appearance. Multiple radiogenomic studies have since found associations between radiomic features and gene expression, highlighting the possible link between biology and radiology [40].

Several studies indicated that skewness is associated with poorer response to radiotherapy and radio-resistance in both lung cancer and other tumor types [34,41,42]. Possible correlations between radiological skewness and tumor biology have been identified in previous studies. Positive skewness was found to be associated with KRAS mutations in NSCLC [43]. Indeed, lung cancers with KRAS mutations are characterized by poor prognosis and treatment resistance and same has been demonstrated in colorectal cancers as well [44,45]. As KRAS is a common mutation in multiple malignancies, the poorer prognosis observed in skewed tumors may be at least partially driven by KRAS mutation [46]. There are multiple mechanisms in which KRAS mutation may lead to poor prognosis, including the induction of a cancer stem cell like phenotype leading to radioresistance [47]. This may be a possible explanation for diminished response towards SBRT in patients with positive skewness in the present study.

Low FOXF2 expression was a newly discovered as an adverse prognostic factor for NSCLC [48]. FOXF2 inhibits proliferation, migration and invasion of tumor cells in vitro and was shown to enhance metastasis [49]. Perez-Morales et al. [39] demonstrated that lung tumor RMS level was correlated with FOXF2 expression. Hence, it is possible that radiological RMS level may reflect underlying FOXF2 expression in the tumor.

Such evidence reinforce the concept that tumor radiomic features might reflect associated underlying biological processes which ultimately affects disease prognosis and radio-sensitivity.

4. Application of machine learning in radiomics

Tumor textural analysis allows intralesional heterogeneity to be represented by myriad of statistical representations. Thus, datasets are typically high dimensional. In the present study, 110 radiomic features were evaluated on 85 analyzed lesions. Hypothesis testing through univariate models would necessitate multiple statistical testing. This leads to higher type I error, hence, Bonferroni adjustment would be required. Furthermore, one particular radiomic feature may not be sufficient for good classification performance. This is evident as the AUC of skewness and RMS were only 0.625 and 0.632, respectively. The optimal decision boundary may lie in higher dimensional space which is better suited for multivariate modelling. SVM allows the decision boundary to lie within a hyperplane in high-dimensional space, as such classification performance would be significantly better [50]. SVM searches for a hyperplane that can maximally separate the outcomes of interest. Multiple studies have demonstrated the ability of SVM to classify high-dimensional clinical dataset for prediction of clinical outcomes such as development of pneumonitis or overall survival following radiotherapy [51,52]. We have demonstrated in this study that SVM was highly accurate in predicting tumor response following SBRT.

Clinical outcomes in this study were unbalanced, as patients with CR (n = 31) and PR (n = 11) were less than NR (n = 43). Despite the unbalance dataset, we were still able to identify two radiomic features that are predictive of NR which are in accordance with previous studies. In datasets with unbalanced class, the accuracy of minority class prediction is often sacrificed by the algorithm to allow higher overall accuracy [53]. Hence while building the SVM model, PR and CR were up-sampled using ADASYN to match the NR class. This approach preserves the statistical distribution of up-sampled class while increasing the statistical power [53]. Similar approach using the Synthetic Minority Oversampling Technique (SMOTE) was demonstrated in previous radiomic studies to be beneficial toward overall classification performance [19].

The application of machine learning models could have great potential in the field of radiomics and precision medicine. Machine learning allows multiple factors to be accounted in a multi-dimensional model which can be useful in the age of precision medicine. The application of radiomics and machine learning in clinical practice requires accurate and reproducible models. As radiological response to SBRT can potentially directly impact survival, the predictive ability of radiomic models can guide decisions to select tumors that are likely responsive to SBRT while opting for alternative therapies for more radioresistant tumors. The optimization of dedicated models could thus allow the introduction of radiomics in the clinical decision-making process.

5. Limitations

1) Here we present a hypothesis generating study aiming to elucidate prognostic radiomic features for future confirmation

Despite being a promising technique to provide additional prognostic information, this study is characterized by several limitations. This retrospective study enrolled a small and heterogeneous cohort which may limit its generalizability. As the study was retrospectively conducted, some factors may not be adequately controlled for. As described previously, the outcomes in each category (CR, PR, and NR) were unbalanced. Although ADASYN can compensate for the imbalance statistically, genuine data points are still more preferable.

Apart from that, we acknowledge that the proportion of patients with NSCLC and non-lung tumors were unbalanced as well (55 vs. 14 patients). However, on a tumor level, there were 62 NSCLC lesions versus 23 non-lung tumors. The level of unbalance was also reduced by ADASYN to a certain extent.

Planning and treatment variability also contribute to limitation of this study. Excessive tumor movement during planning CT may limit the accuracy of feature extraction. This was addressed by recruiting only patients irradiated using ABC technique to reduce tumor movement. It should also be noted that a few fractionation schedules were used within in study which may lead to differential tumor response.

This study had a strict patient selection criteria in which patients with central tumors and patients that cannot tolerate ABC were excluded. This improved the accuracy of the model but at the same time limited its application in “real-world” clinical practice.

We acknowledge that different institutions utilize different imaging protocols and scanners. In addition, use of contrast or injection time may differ across institution. It has been shown that use of contrast may affect certain radiomic features [54]. Non-contrast CT thorax was employed in this study while other studies employed contrast planning CT. This may limit the generalizability and sharing of radiomic models. As such, radiomic models application may be limited to a single institution without the inclusion of multi-institutional data. A universal feature extraction protocol and feature set were also lacking, which hinders inter-study comparison. This can be overcome by data-sharing and cross-validation across different data sets which hopefully will result in more generalizable models in the future.

6. Conclusion

The present study provided a proof of concept, that radiomic features can predict response to SBRT, irrespective of tumor histology, site of primary and initial staging. We also demonstrated that skewness and RMS are potential universal adverse predictive factors, in accordance with multiple previous radiomic studies. It was also shown that there was a trend towards improved overall survival tumors with a radiological response to SBRT. This validates the role of radiological response as a prognostic factor of overall survival.

The clinical implications of the definition of universal radiomic predictive features require additional prospective analysis. Further research is warranted to unravel the biological basis underlying radiomic features and treatment responsiveness.

Multivariate models such as SVM can improve classification performance allowing accurate outcome prediction. Radiomic analysis using such models can potentially be applied for treatment selection. Data sharing and unified feature extraction and processing procedures are eagerly awaited to generate more robust models in the future.