Establishment of a Prediction Model for Overall Survival after Stereotactic Body Radiation Therapy for Primary Non-Small Cell Lung Cancer Using Radiomics Analysis

Simple Summary Lung cancer remains the leading cause of cancer-related mortality worldwide. Although early-stage non-small cell lung cancer (NSCLC) is likely to be controlled with stereotactic body radiation therapy (SBRT), approximately 18% of patients lead to recurrence. The aim of this study was to identify prognostic factors and establish a predictive model for survival outcomes of patients with non-metastatic NSCLC treated with SBRT. Several radiomic features were selected as predictive factors and two prediction models were established from the pre-treatment computed tomography images of 250 patients in the training cohort. One radiomic factor remained a significant prognostic factor of overall survival (OS) (p = 0.044), and one predicting model could estimate OS time (mean: 37.8 months) similar to the real OS time (33.7 months). In this study, we identified one radiomic factor and one prediction model that can be widely used. Abstract Stereotactic body radiation therapy (SBRT) for early-stage non-small cell lung cancer (NSCLC) leads to recurrence in approximately 18% of patients. We aimed to extract the radiomic features, with which we predicted clinical outcomes and to establish predictive models. Patients with primary non-metastatic NSCLC who were treated with SBRT between 2002 and 2022 were retrospectively reviewed. The 358 primary tumors were randomly divided into a training cohort of 250 tumors and a validation cohort of 108 tumors. Clinical features and 744 radiomic features derived from primary tumor delineation on pre-treatment computed tomography were examined as prognostic factors of survival outcomes by univariate and multivariate analyses in the training cohort. Predictive models of survival outcomes were established from the results of the multivariate analysis in the training cohort. The selected radiomic features and prediction models were tested in a validation cohort. We found that one radiomic feature showed a significant difference in overall survival (OS) in the validation cohort (p = 0.044) and one predicting model could estimate OS time (mean: 37.8 months) similar to the real OS time (33.7 months). In this study, we identified one radiomic factor and one prediction model that can be widely used.

Appropriate prediction of recurrence risk before treatment is essential for personalized treatment. If we can select patients with a high possibility of recurrence in advance, we This study was reviewed and approved by the institutional review board and ethics committee (examination number 3372). We retrospectively reviewed an institutional database to identify patients with primary non-metastatic NSCLC who were treated with SBRT with curative intent between April 2002 and March 2022 with or without salvage treatment. Patients who had plain chest CT scans within two weeks prior to SBRT were included. Patients who had never been examined to confirm the clinical outcomes after SBRT were excluded. We detected 358 primary tumors in 338 patients. A total of 250 primary tumors were randomly assigned to the training cohort, and 108 tumors were assigned to the validation cohort.

Radiotherapy
All patients with no pathological diagnosis by biopsy required a continuous increase in the primary tumor size or abnormal accumulation of fluorodeoxyglucose in the tumor before SBRT. Tumor markers (neuron-specific enolase, NSE, and pro-gastrin-releasing peptide, ProGRP) needed to be checked to exclude small cell carcinoma. The patient's forced expiratory volume in one second needed to be more than 750 cc to receive SBRT. In addition, indications for SBRT must be decided by an institutional cancer board consisting of pulmonologists, thoracic surgeons, radiologists, and radiation oncologists.
Six or 10 MV photon linear accelerators were used for SBRT. The radiation technique used was three-dimensional conformal radiotherapy (3D-CRT) or volumetric modulated arc therapy (VMAT). The fixed multi-port technique consisting of 7-12 beams was used as 3D-CRT. A total dose of 42-64 Gy in 4-10 fractions was administered. When the α/β ratio is assumed to be 10, BED was 75-166 Gy. The total dose and fractions were determined by radiation oncologists by considering the location of the tumor and the expected dose of organs at risk. Gross tumor volume (GTV) was contoured on the planning CT images displayed at a lung window level. The clinical target volume (CTV) was defined to be the same as the GTV. An internal margin was added to the CTV to cover the respiratory motion of the tumor and generate the internal target volume (ITV) using four-dimensional CT images. A tracking system was not used. Then, a 5 mm margin in every direction was added to the ITV to generate the planning target volume (PTV). Radiation dose was prescribed to 95% of the PTV, and radiation was administered every weekday.

Clinical Endpoints
As background features, age, Karnofsky performance status (KPS), sex, total radiation dose, radiotherapy (RT) technique (fixed multi-ports, VMAT, or both), primary tumor location, maximum tumor diameter, and SUV-max in PET were investigated.
Overall survival (OS), local relapse-free survival (LRFS), and progression-free survival (PFS) were used as the endpoints. These were counted from the date of SBRT initiation. The patients were usually followed-up and underwent a CT scan every 2-3 months for the first year and repeated 3-6 months thereafter. Local recurrence, which meant a recurrence in the radiation field, was usually identified on PET when SUV-max of the tumor was higher than 2.5 or by biopsy in case of obvious enlargement of the local tumor on the follow-up CT scan. However, some relapses were judged only by continuous enlargement on follow-up CT.

CT Image and Tumor Contouring
All patients underwent plain chest CT for treatment planning. We used Aquilion LB (Canon Medical Systems Corporation, Otawara, Japan) as the CT scanner. The CT tube voltage was 120 kV, the tube current was 350 mA, and the slice thickness was 2 mm. Every primary tumor was manually delineated as the GTV for this study on the expiratory phase of the planning CT displayed with a lung window level using Monaco (Elekta AB, Stockholm, Sweden) by a single expert radiation oncologist. Tumor outline was contoured, and pleura, vessels, heart, bones, and chest wall were excluded from the GTV.

Radiomic Feature Extraction
We used PyRadiomics v3.0.1 (Boston, MA, USA) to extract radiomic features from contoured regions of interest (ROIs) of GTVs [34]. A fixed bin count of 64 was used for discretization of the image gray level. All CT voxels were resampled to 1 × 1 × 1 mm 3 using a B-Spline interpolation function. A total of 107 quantitative features were automatically extracted, including 19 first-order statistics features (intensity histogram, IH), 26 shape-based histogram features, and texture features (gray-level co-occurrence matrix, GLCM, 24 features; gray-level run-length matrix, GLRLM, 16 features; gray-level size-zone matrix, GLSZM, 16 features; neighboring gray-tone difference matrix, NGTDM, 5 features; and gray-level dependence matrix, GLDM, 14 features). All quantitative features are listed in Table S1. From these original radiomic features, 744 quantitative features were obtained using wavelet transformation. Coiflets 1, which is one of the methods of wavelet transformation, was used in x, y and z axes in this study. As a result, 8 wavelettransformed images (HHH, HHL, HLH, HLL, LHH, LHL, LLH, and LLL) were obtained, where H represents high frequency and L represents low frequency. Mathematical definitions of these radiomic features have previously been described and are available at https://pyradiomics.readthedocs.io/en/latest/features.html (accessed on 11 June 2022).

Statistical Analysis
The Kaplan-Meier method was used for survival time analysis. Clinical features of continuous variables consisted of maximum tumor diameter, age, KPS, SUV-max in PET, and total radiation dose, and radiomic features were divided into two groups according to the median values in the training cohort. Univariate analysis was performed using the log-rank test to assess the difference in clinical outcomes between the high-and low-value groups. Multivariate analysis was conducted using the Cox proportional hazards model among the factors with a significant difference in univariate analysis. The stepwise method with the Bayesian information criterion was used for variable selection in the training cohort (direction = "backward/forward"). We created a prediction score from the results of the Cox proportional hazards model and multiple linear regression analysis of the training cohort. The selected radiomic factors and prediction scores were validated in the validation cohort. All statistical analyses were performed using R v4.1.3. The results of the univariate analysis of radiomic features were considered statistically significant at p < 0.01, whereas the other results were considered statistically significant at p < 0.05.

Patient Characteristics
Patients were randomly divided into training (n = 250) and validation (n = 108) groups. As shown in Table 1, there was no bias in the proportions or means of both the patient and tumor background factors in either group. In the training group, 68% were male, mean age was 77 (range, 42-93) years, 84% had peripheral lesions, mean maximum tumor diameter was 20.5 mm (4.0-90.0), mean total dose was 52.5 Gy (42-64), 70% received VMAT, and 12.8% received salvage treatment after recurrence was recognized. No patient underwent adjuvant treatment after SBRT. The median follow-up time for patients who survived at the last follow-up of the training cohort and the validation cohort was 33.0 (0.9-167.9) and 33.9 (0.5-149.4) months, respectively.

Univariate Analysis
As shown in Table 2, as background factors affecting OS, KPS, maximum tumor diameter, sex, peripheral lesions, and salvage treatment showed a significant difference of p < 0.05 in the univariate analysis using the log-rank test. Regarding LRFS, KPS, maximum tumor diameter, sex, peripheral lesions, consolidation/tumor (C/T) ratio, salvage treatment, and SUV-max showed a significant difference of p < 0.05 in the univariate analysis. For PFS, KPS, peripheral lesions, and salvage treatment showed a significant difference of p < 0.05 in the univariate analysis.  Among the 744 radiomic features, 70 had a significant difference of p < 0.01 in OS between the high-and low-value groups in the training cohort in the log-rank test. For LRFS, 223 factors, and for PFS, 210 factors showed a significant difference of p < 0.01.

Prediction Score
From the results of the Cox proportional hazards regression model, we created a prediction score, which was named the Cox-score, calculated by the following formula using the coef value as a coefficient.
Cox For the [x] value, "+1" is substituted if the value of x is higher than the median value, and "−1" is substituted if the value is lower than the median value.
We estimated survival time using multiple linear regression analysis in the training cohort. The estimated survival time was calculated using the following formula.
Estimated For the [x] value, "+1" is substituted if the value of x is higher than the median value, and "0" is substituted if the value is lower than the median value.

Validation
As a result of the log-rank test of 108 cases of the validation cohort for each of the radiomic factors remaining in the multivariate analysis, only one factor of "LargeAreaEm-phasis_LHH" showed a significant difference in OS (p = 0.044, 5-year OS of 70.7% vs. 50.3%) (Figure 1). "LargeAreaEmphasis_LHH" was a parameter that was categorized in GLSZM features. Regarding LRFS and PFS, none of the factors that showed a significant difference in the multivariate analysis of the training cohort showed a significant difference in the univariate analysis of the validation cohort.  Regarding the Cox-score, the value for each case was calculated using the distribution of the radiomic features of 108 cases in the validation group. As a result, 38 cases with a score of less than −0.1 were divided into a low value group, 36 cases with a score of +1.0 or more were divided into a high value group, and the other 34 cases were divided into a median value group for OS. In the log-rank test, the difference in OS between the three groups was not significant (p = 0.63, 5-year OS of 62.6% in the high score group, 59.6% in the median, and 65.4% in the low score group). The difference in LRFS between the two groups, if +1.0 was the cutoff value, was not significant (p = 0.086, 5-year OS of 60.8% in the high, and 82.8% in the low). The difference in PFS between the two groups, if −1.0 was the cutoff value, was not significant (p = 0.070, 5-year OS of 52.0% in the high, and 94.7% in the low).
We compared the estimated survival time calculated using the formula from multiple linear regression analysis with the real survival time in the validation cohort. As a result of the estimated OS time, the difference in the average values was 4.1 months (33.7 months in the real survival time vs. 37.8 months in the estimated survival time) and the t value and p value were 1.619 and 0.11 by the paired t-test, respectively (Figure 2). Regarding LRFS and PFS, the p value was less than 0.0001 using the paired t-test. Figure 1. Impact of the radiomic factor named "LargeAreaEmphasis_LHH" on overall survival in the validation cohort. Patients were divided into two groups according to the median value of "LargeAreaEmphasis_LHH" (low value group and high value group). The univariate analysis using the log-rank test was conducted (p = 0.044, 5-year OS of 70.7% vs. 50.3%).
Regarding the Cox-score, the value for each case was calculated using the distribution of the radiomic features of 108 cases in the validation group. As a result, 38 cases with a score of less than −0.1 were divided into a low value group, 36 cases with a score of +1.0 or more were divided into a high value group, and the other 34 cases were divided into a median value group for OS. In the log-rank test, the difference in OS between the three groups was not significant (p = 0.63, 5-year OS of 62.6% in the high score group, 59.6% in the median, and 65.4% in the low score group). The difference in LRFS between the two groups, if +1.0 was the cutoff value, was not significant (p = 0.086, 5-year OS of 60.8% in the high, and 82.8% in the low). The difference in PFS between the two groups, if −1.0 was the cutoff value, was not significant (p = 0.070, 5-year OS of 52.0% in the high, and 94.7% in the low).
We compared the estimated survival time calculated using the formula from multiple linear regression analysis with the real survival time in the validation cohort. As a result of the estimated OS time, the difference in the average values was 4.1 months (33.7 months in the real survival time vs. 37.8 months in the estimated survival time) and the t value and p value were 1.619 and 0.11 by the paired t-test, respectively (Figure 2). Regarding LRFS and PFS, the p value was less than 0.0001 using the paired t-test.

Discussion
In this study, we extracted radiomic features from the GTV on SBRT planning CT of patients with non-metastatic primary NSCLC and selected features correlated with survival outcomes. We established a prediction model using these selected features. The selected features and prediction models were validated in an independent external cohort. "LargeAreaEmphasis_LHH" categorized in GLSZM features remained as a factor significantly correlated with OS in the validation cohort. The estimated OS time calculated using the formula from the multiple linear regression analysis was similar to the actual OS time in the validation cohort. There was no correlation between the radiomic features and LRFS or PFS.
Prediction of recurrence risk before treatment for early-stage NSCLC is so important. The treatment intensity such as the total radiation dose, fractions, and radiation field can be increased for high-risk patients. Another measure to increase the treatment intensity is combined therapy. Ernani et al. [35] showed the promising OS data of SBRT and adjuvant chemotherapy for patients with tumors less than 4 cm and node-negative NSCLC. Use of immune checkpoint inhibitors as an adjuvant therapy after curative chemoradiotherapy for advanced stage NSCLC have rapidly spread worldwide [36,37]. Several studies have also reported the effectiveness of immune checkpoint inhibitors together with SBRT for high-risk patients with early-stage NSCLC [38]. The selected radiomic features and OS prediction model derive in the present study may contribute to selecting patients who benefit from these treatment options.
Several studies have investigated the prognostic factors of patients who received SBRT for early-stage NSCLC using radiomics analysis. Oikonomou et al. [32] examined CT and PET-derived radiomic features in 150 patients with non-metastatic NSCLC treated with SBRT. They examined several principal components, consisting of six to eight radiomic features. Four principal components were found to be significantly associated with the survival outcomes. In the combined analysis with SUV-max, radiomic features remained the only predictors of OS, disease specific survival, and regional control. Huynh et al. [30] analyzed 113 patients with stage I-II NSCLC who were treated with SBRT. Among the 1605 radiomic features extracted from the tumors in the free breathing CT images, one radiomic feature was a significant prognostic factor for distant metastasis and four radiomic features were prognostic for OS. Starkov et al. [33] retrospectively analyzed 116 patients treated with SBRT for biopsy-confirmed primary NSCLC at a single institution. They showed that quantitative image features from NSCLC nodules on pre-

Discussion
In this study, we extracted radiomic features from the GTV on SBRT planning CT of patients with non-metastatic primary NSCLC and selected features correlated with survival outcomes. We established a prediction model using these selected features. The selected features and prediction models were validated in an independent external cohort. "LargeAreaEmphasis_LHH" categorized in GLSZM features remained as a factor significantly correlated with OS in the validation cohort. The estimated OS time calculated using the formula from the multiple linear regression analysis was similar to the actual OS time in the validation cohort. There was no correlation between the radiomic features and LRFS or PFS.
Prediction of recurrence risk before treatment for early-stage NSCLC is so important. The treatment intensity such as the total radiation dose, fractions, and radiation field can be increased for high-risk patients. Another measure to increase the treatment intensity is combined therapy. Ernani et al. [35] showed the promising OS data of SBRT and adjuvant chemotherapy for patients with tumors less than 4 cm and node-negative NSCLC. Use of immune checkpoint inhibitors as an adjuvant therapy after curative chemoradiotherapy for advanced stage NSCLC have rapidly spread worldwide [36,37]. Several studies have also reported the effectiveness of immune checkpoint inhibitors together with SBRT for high-risk patients with early-stage NSCLC [38]. The selected radiomic features and OS prediction model derive in the present study may contribute to selecting patients who benefit from these treatment options.
Several studies have investigated the prognostic factors of patients who received SBRT for early-stage NSCLC using radiomics analysis. Oikonomou et al. [32] examined CT and PET-derived radiomic features in 150 patients with non-metastatic NSCLC treated with SBRT. They examined several principal components, consisting of six to eight radiomic features. Four principal components were found to be significantly associated with the survival outcomes. In the combined analysis with SUV-max, radiomic features remained the only predictors of OS, disease specific survival, and regional control. Huynh et al. [30] analyzed 113 patients with stage I-II NSCLC who were treated with SBRT. Among the 1605 radiomic features extracted from the tumors in the free breathing CT images, one radiomic feature was a significant prognostic factor for distant metastasis and four radiomic features were prognostic for OS. Starkov et al. [33] retrospectively analyzed 116 patients treated with SBRT for biopsy-confirmed primary NSCLC at a single institution. They showed that quantitative image features from NSCLC nodules on pre-treatment CT images were correlated with OS after SBRT. Similar to the present study, no significant relationship was observed between radiomic features and survival outcomes except for OS. No external validation cohort was established in any of these studies. In contrast, Dissaux et al. [28] built a prediction model of local control derived from PET among a training set of 64 earlystage NSCLC patients treated with SBRT. The model maintained a significant correlation with local control in the validation set of 23 patients. The fact that patients in the training set were collected from two institutions and those in the validation set were collected from another single institution may have affected the results of the study. Yu et al. [39] retrospectively extracted radiomic features from contrast-enhanced CT images of a training cohort (n = 147) and an independent validation cohort (n = 295) of patients with stage I NSCLC. Although the cohort size was large, there was an important difference between the two cohorts. The prediction model consisting of two radiomic features was significantly associated with both OS and distant metastasis in the validation cohort. Patients in the training cohort were treated with surgery, and those in the validation cohort were treated with SBRT. In addition, the patients in the validation cohort were significantly older than those in the training cohort, which seemed to be due to the difference in the indication for treatment. These differences may have affected the accuracy of the validation. In contrast, the training and validation cohorts were randomly divided in the present study. Thus, the clinical background features were not significantly different between the two cohorts. Franceschini et al. [29] also analyzed radiomic features in the treatment planning CT images of 70 patients treated with SBRT for early-stage NSCLC for training and 32 patients for validation. The two cohorts were randomly divided as in the present study. A set of 45 textural features was extracted from the tumor volumes on the treatment planning CT images. They built a predictive model that was prognostic of PFS and disease-specific survival. Although these studies included a validation process, the sample size may not have been large enough. In our study, we examined much more NSCLC nodules, which were randomly divided into training and validation cohorts to minimize the risk of bias. In addition to background factors, crude OS, LRFS, and PFS were similar between the two cohorts. Several radiomic features remained as factors correlated with PFS and LRFS in the multivariate analysis of the training cohort; however, none of them were found to be prognostic factors in the external validation process. Unless validation is conducted, these radiomic features could be identified as prognostic factors.
The present study has some limitations. First, although PET or biopsy was recommended to identify local recurrence, there were some cases in which local recurrence was judged using only CT. Thus, there may have been an overdiagnosis of local recurrence. Second, all patients were recruited from a single institution and all tumors were contoured by a single radiation oncologist. External validation at various institutions is expected to be conducted to develop the established prediction model into a widely usable model. Third, the fact that all GTVs were manually contoured may have reduced reproducibility of this study. Homogeneity of this study was kept since all contouring tasks were conducted by a single expert radiation oncologist.

Conclusions
We found that the radiomic features categorized as GLSZM features derived from the GTV on the pre-treatment CT image were an independent prognostic factor of OS in non-metastatic NSCLC patients treated with curative SBRT. Our findings suggest that the estimated OS time derived from the multiple linear regression analysis was similar to the real OS time. These findings may contribute to selecting patients who benefit from increasing the treatment intensity by, for example, increasing radiation dose and using combined therapy.
The selected radiomic factor and OS prediction model derived in this study can be widely used in the future. External validation at various institutions is expected.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/cancers14163859/s1, Table S1: List of all 107 radiomic features. Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.