Prognostic Impact of Longitudinal Monitoring of Radiomic Features in Patients with Advanced Non-Small Cell Lung Cancer

Tumor growth dynamics vary substantially in non-small cell lung cancer (NSCLC). We aimed to develop biomarkers reflecting longitudinal change of radiomic features in NSCLC and evaluate their prognostic power. Fifty-three patients with advanced NSCLC were included. Three primary variables reflecting patterns of longitudinal change were extracted: area under the curve of longitudinal change (AUC1), beta value reflecting slope over time, and AUC2, a value obtained by considering the slope and area over the longitudinal change of features. We constructed models for predicting survival with multivariate cox regression, and identified the performance of these models. AUC2 exhibited an excellent correlation between patterns of longitudinal volume change and a significant difference in overall survival time. Multivariate regression analysis based on cut-off values of radiomic features extracted from baseline CT and AUC2 showed that kurtosis of positive pixel values and surface area from baseline CT, AUC2 of density, skewness of positive pixel values, and entropy at inner portion were associated with overall survival. For the prediction model, the areas under the receiver operating characteristic curve (AUROC) were 0.948 and 0.862 at 1 and 3 years of follow-up, respectively. Longitudinal change of radiomic tumor features may serve as prognostic biomarkers in patients with advanced NSCLC.

Quantitative imaging features may provide comprehensive information on tumor phenotypes and microenvironement 9,10 . Changes in quantitative features are considered a biomarker to accurately assess treatment response and clinical outcomes 11 . Recently, delta-radiomics, a time dependent metric comprised of quantitative features extracted from medical images acquired during the course of treatment has been considered as an emerging biomarker for assessing therapeutic responses and clinical outcomes 12 . We hypothesized that variables capable of reflecting various longitudinal changes in quantitative image features extracted from serial imaging over the course of treatment would be capable of superior prognosis prediction in NSCLC patients. The purpose of this study was to develop biomarkers based on longitudinal changes in quantitative radiomic features in patients with advanced NSCLC, and assess the prognostic power of these biomarkers.
Development of radiologic variables reflecting longitudinal volume change. AUC1, AUC2, and beta value (Fig. 1) were calculated using the volume change of serial CTs from baseline to PD. The beta value exhibited a good correlation with the six patterns of volume change (r = 0.834, p < 0.000). Furthermore, AUC2 reflecting under the area and slope during the follow-up period was well correlated with the six patterns of volume change (r = 0.848, p < 0.000) as well as beta value (r = 0.752, p < 0.000, Fig. 2). Therefore, AUC2 was considered to be a marker that reflects changes during the follow-up in patients with advanced NSCLC and was used for longitudinal analyses of the 23 quantitative CT features, including volume.
For the prediction of overall survival time, we calculated the cutoff of AUC1, AUC2, and beta value based on volume. When we selected cutoff values of 1.601 for AUC1, −0.738 for AUC2, and −0.176 for the beta value, the www.nature.com/scientificreports www.nature.com/scientificreports/ AUC2 and beta values showed significant differences in overall survival time (p = 0.035 for AUC2, and p = 0.029 for beta value (Fig. 3). However, AUC1 did not exhibit a statistically significant difference with regards to overall survival (p = 0.086).
prognostic prediction model using longitudinal change of radiomic features. We extracted 23 quantitative features including volume from baseline and follow-up CTs, and calculated the AUC2 of the quantitative CT features extracted from baseline to PD. Based on the cut-off values of quantitative features on baseline CT and AUC2 of quantitative values, all quantitative features were divided into two groups and univariate and multivariate cox regression analyses were performed to identify risk factors to predict overall survival time. Eighteen features were selected in the univariate analysis (supplementary Table 1). In our study, age, sex, the Eastern Cooperative Oncology group (ECOG) performance status, smoking status, type of EGFR mutation, CNS metastasis, and local treatment for CNS metastasis were not associated with overall survival time (p > 0.05). Six of the 18 features were selected by stepwise selection. Kurtosis of positive pixel values (p = < 0.000, HR, 7.992; 95% CI, 2.759-23.148) and surface area (p = 0.001, HR, 0.192; 95% CI, 0.075-0.493) on baseline CT, and AUC2 of density (p = < 0.000, HR, 0.067; 95% CI, 0.021-0.220), skewness of positive pixel values (p = 0.010, HR, 5.142; 95% CI, 1.490-17.747) and entropy at inner portion (p = 0.001, HR, 6.196; 95% CI, 2.086-18.406) were associated with overall survival in multivariate cox analysis in patients with advanced NSCLC ( Table 2).
The prediction model of overall survival was defined as a linear combination of the five selected features and their regression coefficients. For the prediction model, we constructed a time-dependent ROC curve from which the AUROC was calculated. At 1 and 3 years of follow-up, the AUROCs for overall survival were 0.948, and 0.862 (Fig. 4).  The six patterns of volume changes were follows: (1) reduction only but progressive disease due to nontarget lesions, (2) slow progression after rapid response, (3) rapid progression after rapid response, (4) slow progression after slight reduction, (5) rapid progression after slight reduction, and (6) sequential progression.

Discussion
Radiomics studies of NSCLC are usually obtained from images at a single time point. Our studies extracted quantitative features from serial images obtained during the course of treatment in patients with advanced NSCLC. As a result of analyzing longitudinal changes in volume, we identified the variable AUC2 that reflects the slope and area of longitudinal change. Specifically, AUC2 of volume was helpful in predicting overall survival time. Based on this result, we also obtained the AUC2 for 23 quantitative features from serial CT images. Kurtosis of positive pixel values and surface area from baseline CT, AUC2 of density, skewness of positive pixel values, and entropy at inner portion were all associated with overall survival.   www.nature.com/scientificreports www.nature.com/scientificreports/ Imaging plays a central role in the evaluation of therapeutic response. Decrease in size, vascularity and metabolism, and cavitation within the tumor are well known typical responses to treatment 13 . However, although tumor size can increase with intralesional and/or perilesional hemorrhage, atypical responses consisting of decreased attenuation and metabolism of tumor may also be seen as the tumor becomes homogenous 13 .
As the use of targeted molecular therapy and immunotherapy increases in the era of personalized cancer treatment, morphological changes may be more important than changes in tumor size, and increases in the size of new lesions before dramatic shrinkage may also occur in relation to the immune response 14 . Patients treated with EGFR TKIs can experience a dramatic initial shrinkage of tumor followed by slow progression over several months, a finding that suggests the tumor cells remain sensitive to EGFR TKIs 2,14 . A recent study reported that 88% of patients continued therapy beyond progression according RECIST 15 . Evaluating therapeutic response as a more continuous variable using DepoR has also begun to be explored, and DepoR is used not only in lung cancer, but in other cancers as well including colorectal cancer, breast cancer, and multiple myeloma 16 . A recent study showed that a greater DepoR is associated with better clinical outcome in NSCLC patients receiving an anaplastic lymphoma kinase (ALK) inhibitor or a programed cell death-1 (PD-1) inhibitor 5 . On the other hand, another study reported that DepoR does not correlate with clinical outcome 8 . Therefore, an approach based on a single marker is not sensitive to assessment of therapeutic response due to phenotypic and genomic heterogeneity between lesions, even within an individual tumor 17,18 .
Radiomics is defined as quantitative mapping that converts medial images to automatically extracted quantitative features and provides information related to the comprehensive pathophysiology of a tumor 12,19 . The use of radiomics in lung cancer studies is focused on patient outcomes, nodule discrimination, and assessment of genetic mutations 20 . Whereas most radiomics studies are performed using pretreatment images, the concept of delta-radiomics was introduced to analyze quantitative features from images acquired over the course of treatment by adding a time component to existing radiomic features 12 . Delta-radiomics may be used as indicator to predict response in colorectal liver metastasis, metastatic lung cancer, esophageal cancer, and NSCLC [21][22][23][24][25] . Previous studies have shown that tumor volume is a key predictor reflecting histological tumor aggressiveness and poor prognosis in patients with NSCLC 26,27 . Volumetric measurements that are able to capture tumor changes along the z-axis may be useful for more sensitively and specifically monitoring disease change over time [28][29][30] . Therefore, in our study, we examined longitudinal changes in tumor volume. We measured additional variables reflecting longitudinal change, and from those variable calculated corresponding AUC1 values reflecting the area of volume change, beta value reflecting slope of change, and AUC2 reflecting area and slope. AUC2 exhibited a good correlation with curve pattern (r = 0.848, p < 0.000) and beta value (r = 0.752, p < 0.000), whereas the AUC2 of tumor volume was associated with overall survival (p = 0.035). Thus, AUC2 may be a new biomarker that reflects longitudinal changes in lesions and predicts overall survival time and treatment response.
Size and volume of a tumor do not reflect changes in tumor heterogeneity and genetic profiles 23 . However, radiomic features can provide more comprehensive information about tumor makeup. A study showed that a novel set of quantitative image features, based on heterogeneities of tumor physiology, was helpful for early prediction of treatment outcome 31 . Therefore, we obtained AUC2 values for 23 quantitative features, including volume, that were previously reported to be associated with lung cancer [32][33][34][35][36][37] . We examined the association between these 23 pretreatment quantitative features and 23 AUC2s obtained from serial images with overall survival. The result of our study showed that kurtosis of positive pixel values (p = <0.000) and surface area (p = 0.001) on baseline CT, and AUC2 of density (p < 0.000), skewness of positive pixel values (p = 0.010), and entropy at inner portion (p = 0.001) were all associated with overall survival according to a multivariate analysis in patients with advanced NSCLC.
Targeted therapies developed in recent years can cause tumor necrosis, hemorrhage, or cavitation, and radiomic features based on histograms may be useful for predicting these pathologic response 27 . In patients with a gastrointestinal stromal tumor, treatment with imatinib results in a decrease in tumor density such that the appearance becomes homogenous on CT, even though the decrease in tumor size is minimal at an early posttreatment time point 38 . In our study with advanced NSCLC, density change was associated with increased survival, and it has been suggested that increased homogeneity during the course of treatment is related with therapeutic response and prognosis. Two studies reported that a reduction of entropy reflecting tumor heterogeneity as a parameter to evaluate treatment response in colorectal liver metastases and renal cell carcinoma 24,25 . In our study, we found that the AUC2 of entropy at inner portion was successful for predicting overall survival time in advanced NSCLC. In addition, the AUC2 of the skewness of positive pixel value was found to be predictor of overall survival. Skewness is the distribution pattern of CT attenuation values, and decreased skewness reflects decreased enhancement due to decreased neovascularization during the course of treatment 34 . Therefore, based on our results, changes in tumor heterogeneity were related to overall survival time.
Morphological features reflecting the physical characteristics of tumors are important features related to aggressiveness 18,39 . Irregularity of tumor shape reflecting non-uniform growth of tumors is associated with worse overall survival in patients with lung squamous cell carcinoma and lung adenocarcinoma 32,40 . In our study, surface area on baseline CT was associated with overall survival time in patients with advanced lung cancer, and morphologic features of tumor were considered important predictors related to clinical outcomes.
This study had several limitations. First, the sample size of the present study was small, and our study included patients in various treatment course. However, inclusion of patients with various treatment course may provide a better, more comprehensive assessment of overall survival of patients with advanced lung cancer. A second limitation of this study was that we did not externally validate our results. Therefore, subsequent studies with a large population and external validation are necessary for validation of our results. Third, tumor volume and radiomic features are susceptible to respiration due to collapse and stretching of the tumor and peritumoral lung parenchyma 41,42 . In order to minimize this effect, we excluded patients who had severe respiratory difference in chest CTs or had concurrent disease such as effusion. Fourth, we selected target lesions according to RECIST (version (2019) 9:8730 | https://doi.org/10.1038/s41598-019-45117-y www.nature.com/scientificreports www.nature.com/scientificreports/ 1.1) criteria. However, current RECIST guidelines may not reflect the overall body tumor load in patients with metastatic cancer. In addition, selection of target lesions is crucial and the result of response classification could be changed by target lesion selection 43,44 . The final limitation of our study was that the ROIs were drawn manually by a single operator. As tumor segmentation is the most critical process, and it is well recognized that interoperator variability of manual segmentation of lesion is high 9 .

Conclusion
Our study showed that AUC2, a value that reflects the slope and area of longitudinal change calculated from radiomic measurements, is associated with patient survival. Quantitative analysis is a noninvasive technique reflecting comprehensive information such as tumor heterogeneity and genetic profile. In addition, changes in these radiomics features during serial follow-up may provide comprehensive information about response to treatment and prognosis. Therefore, the results of our study suggest that longitudinal change of radiomic tumor features may be useful as prognostic biomarkers in patients with advanced NSCLC.

Materials and Methods
study population. From January 2012 to October 2014, we recruited 80 patients with advanced or recurrent NSCLC with confirmed EGFR mutation treated with a tyrosine kinase inhibitor (TKI) and who had been evaluated for therapeutic response on computed tomography (CT) until PD on treatment. Treatment response was assessed according to Response Evaluation Criteria in Solid Tumors (RECIST version 1.1). There was no change in TKI type during the course of treatment. Twenty-seven patients were excluded from our study based on the following exclusion criteria: (1) patients who underwent surgical resection (n = 9), (2) patients without stage IV disease (n = 1), and (3) patient without available serial CTs (n = 17). Thus, a total of 53 patients with advanced NSCLC were included in this retrospective study. The study was approved by the Institutional Review Board (IRB) of Samsung Medical Center (IRB number 2015-10-108), and the requirement for informed consent was waived.
Image acquisition. All helical CT images were obtained with a 64 detector-row (LightSpeed VCT; GE Healthcare, Waukesha, WI, USA) CT scanner using the following parameters: detector collimation, 1.25 or 0.625 mm; field of view, 36 cm; 125 mA; 120 kVp; beam width, 10-20 mm; beam pitch, 1.375-1.5; section thickness 2.5 mm; and matrix, 512 × 512 mm. All patients underwent chest CT at full inspiration through breath hold to minimize the effect of the tumor motion due to breathing. Chest CT scanning was obtained 90 seconds after the administration of contrast material. A total of 1.5 mL/kg (body weight) Iomeron 300 (Iomeprol, 300 mg iodine/mL; Bracco; Milan, Italy) was injected at an infusion rate of 3 mL/s using a power injector (MCT Plus; Medrad; Pittsburgh, PA, USA). Image data were reconstructed with a soft-tissue algorithm for mediastinal window ranges and a bone algorithm for lung window images. Both mediastinal (width, 300 Hounsfield units [HU]; level, 20 HU) and lung (width, 1500 HU; level −700 HU) window images were displayed for tumor assessment. Chest CT images were obtained every two cycles (8 weeks) during the course of treatment.  Table S2). In all available CTs, volume referred to the measurement of the sum of the target lesions and other quantitative features were measured from the volume weighted average of the target lesions.
We extracted three variables reflecting various patterns of longitudinal change of radiomic features including volume, namely, area under the curve (AUC1), beta value, and AUC2 (Fig. 1). AUC1 is a variable that represents the area of the longitudinal change of a quantitative feature. The beta value refers to the slope calculated by linear regression over time, that is, overall slope of change during follow-up. Third, AUC2 is a value obtained by considering the slope and area of the longitudinal changes of quantitative features. We investigated the relationship between the three variables obtained from the volume of target lesions and six patterns visually divided based on the patterns of longitudinal volume change. The six patterns reflecting longitudinal change of volumes were as follows: 1) reduction only but PD due to non-target lesions, 2) slow progression after rapid response, 3) rapid progression after rapid response, 4) slow progression after slight reduction, 5) rapid progression after slight reduction, 6) and sequential progression. statistical analysis. Analyses were performed using SAS version 9.4 (SAS institute, Cary, NC) and R 3.3.1 (Vienna, Austria; http://www.R-project.org/). A p value less than 0.05 was considered statistically significant.
Correlation between novel markers and pattern of longitudinal volume change. AUC1, beta value and AUC2 were obtained by the following methods. The regression coefficient of radiomic feature's value over time was estimated by linear regression and the estimated regression coefficient is denoted as the beta value. The area of radiomic feature from the base line to the next time point was calculated using the formula for the area of a trapezoid. In this way, the area was calculated sequentially up to PD time point. AUC1 represents the sum of all calculated the areas. If the difference of radiomic feature's values between the current time point and the next time point are positive, add it's area value, otherwise subtract it's area value. The total area was calculated by this method, which is denoted as AUC2. The correlations between AUC1, beta value, AUC2 and the six patterns obtained from the (2019) 9:8730 | https://doi.org/10.1038/s41598-019-45117-y www.nature.com/scientificreports www.nature.com/scientificreports/ longitudinal change of volume were calculated by Spearman's rank correlation coefficient. The cut-off values that best predicted overall survival time based on AUC1, beta value, and AUC2 of volume were selected as the point with the most significant log-rank p-value for all possible cut-off points. Overall survival curves of AUC1, AUC2, and beta value based on cut-off values were estimated using the Kaplan-Meier method.
Prognostic prediction model using multivariate analysis. We calculated AUC2 for each of the 23 quantitative features. The cut-off values that best predicted overall survival time for the 23 baseline quantitative features and AUC2 of 23 quantitative features were selected. Univariate and multivariate analyses of the 23 baseline quantitative features with age, sex, ECOG performance state, smoking status, type of EGFR mutation, and AUC2 of the 23 quantitative features were performed using Cox's proportional hazards model to identify risk factors associated with overall survival. Univariate analysis was used to select 18 features with a p value less than 0.05. Spearman's correlation analysis of 18 features and longitudinal data represented by AUC2 was performed. Some features have a positive or negative correlation. We looked at the variance inflation factor (VIF) to see the collinearity between the features included in the multivariable cox regression. The VIF of all variables was less than and there was no collinearity. To identify significant features on Cox's proportional hazard model, we considered stepwise selection with Akaike's Information Criterion (AIC). For significant features of overall survival, a time-dependent receiver operating characteristic (ROC) curve was constructed, and the area under the ROC curve (AUROC) was calculated.