Peritumoral radiomics features predict distant metastasis in locally advanced NSCLC

Purpose Radiomics provides quantitative tissue heterogeneity profiling and is an exciting approach to developing imaging biomarkers in the context of precision medicine. Normal-appearing parenchymal tissues surrounding primary tumors can harbor microscopic disease that leads to increased risk of distant metastasis (DM). This study assesses whether computed-tomography (CT) imaging features of such peritumoral tissues can predict DM in locally advanced non-small cell lung cancer (NSCLC). Material and methods 200 NSCLC patients of histological adenocarcinoma were included in this study. The investigated lung tissues were tumor rim, defined to be 3mm of tumor and parenchymal tissue on either side of the tumor border and the exterior region extended from 3 to 9mm outside of the tumor. Fifteen stable radiomic features were extracted and evaluated from each of these regions on pre-treatment CT images. For comparison, features from expert-delineated tumor contours were similarly prepared. The patient cohort was separated into training and validation datasets for prognostic power evaluation. Both univariable and multivariable analyses were performed for each region using concordance index (CI). Results Univariable analysis reveals that six out of fifteen tumor rim features were significantly prognostic of DM (p-value < 0.05), as were ten features from the visible tumor, and only one of the exterior features was. Multivariablely, a rim radiomic signature achieved the highest prognostic performance in the independent validation sub-cohort (CI = 0.64, p-value = 2.4×10−5) significantly over a multivariable clinical model (CI = 0.53), a visible tumor radiomics model (CI = 0.59), or an exterior tissue model (CI = 0.55). Furthermore, patient stratification by the combined rim signature and clinical predictor led to a significant improvement on the clinical predictor alone and also outperformed stratification using the combined tumor signature and clinical predictor. Conclusions We identified peritumoral rim radiomic features significantly associated with DM. This study demonstrated that peritumoral imaging characteristics may provide additional valuable information over the visible tumor features for patient risk stratification due to cancer metastasis.


Material and methods
200 NSCLC patients of histological adenocarcinoma were included in this study. The investigated lung tissues were tumor rim, defined to be 3mm of tumor and parenchymal tissue on either side of the tumor border and the exterior region extended from 3 to 9mm outside of the tumor. Fifteen stable radiomic features were extracted and evaluated from each of these regions on pre-treatment CT images. For comparison, features from expert-delineated tumor contours were similarly prepared. The patient cohort was separated into training and validation datasets for prognostic power evaluation. Both univariable and multivariable analyses were performed for each region using concordance index (CI).

Results
Univariable analysis reveals that six out of fifteen tumor rim features were significantly prognostic of DM (p-value < 0.05), as were ten features from the visible tumor, and only one of the exterior features was. Multivariablely, a rim radiomic signature achieved the highest prognostic performance in the independent validation sub-cohort (CI = 0.64, p-value = 2.4×10 −5 ) significantly over a multivariable clinical model (CI = 0.53), a visible tumor radiomics model (CI = 0.59), or an exterior tissue model (CI = 0.55). Furthermore, patient stratification by the combined rim signature and clinical predictor led to a significant improvement on the clinical predictor alone and also outperformed stratification using the combined tumor signature and clinical predictor. PLOS

Introduction
Lung cancer remains the leading cause in cancer-related mortality worldwide [1]. Histologically, adenocarcinoma represents the most common type of non-small cell lung cancer (NSCLC). Locally advanced NSCLC patients represent about 30% of newly diagnosed lung cancer [2]. These patients typically receive a combination of surgery, chemotherapy, and radiation therapy [3]. Despite these treatment approaches, the survival rate of these patients is limited to~25% at five years due to disease progression [4,5]. The limitations of the current treatment approach necessitate novel prognosticators that allow for further stratification of different risk groups and more refined therapeutic strategies. Quantitative imaging has been increasingly employed to assess treatment response to cancer therapy. Especially for lung cancers, CT imaging remains the modality of choice as it is noninvasive, renders anatomical details in high resolution and can quickly capture patient thoracic anatomy so that artifacts due to respiratory motion can be minimized. Planning CT images are routinely acquired in lung cancer patients prior to radiation therapy. Recently, numerous studies have shown imaging-based radiomic features can quantify tumor heterogeneity and hold potential for their application as clinical biomarker for patient stratification [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23]. In particular, radiomic studies have shown CT-derived image features may be prognostic for distant metastasis (DM) and treatment responses in NSCLC [13,14,17,24].
Previous studies have predominantly investigated the association between clinical outcomes and radiomic features within the primary tumor volume [24][25][26][27]. However, recent cancer research has shown evidences that extratumoral lung parenchymal tissues surrounding the primary tumor can become involved as cancer infiltrates and metastasizes. Pathological studies have demonstrated that lung tumor can spread through blood and lymphatic vasculature as well as airspaces in lung parenchyma [28][29][30][31][32][33][34][35] and that extratumoral cancerous presence may lead to worse clinical performance. In all aforementioned modes of cancer spreading, study results have consistently found significantly stronger association with distant or local recurrences for the extratumoral cancerous presence than their intratumoral counterparts [28,31,33]. Thus, we hypothesized that tumor metastatic progression may manifest itself in the imaged peritumoral tissue characteristics and the underlying relationship may be explored using radiomics based profiling of the normal-appearing tissue beyond the identified tumor region. Given the lack of biomarkers for DM in NSCLC, an understanding of the peritumoral tissue radiomics as an imaging biomarker may provide additional information to the existing approaches that only quantify characteristics within the visible tumor volume for identifying patients at higher risk and facilitating improved treatment design.
In this study, we present a radiomics investigation on the association of DM with peritumoral tissues in a cohort of 200 adenocarcinoma NSCLC patients. For clinical utilization, their prognostic performances were compared to the tumor-only radiomic features and clinical factors.

Patient characteristics
This study was conducted under Dana-Farber/Harvard Cancer Center IRB protocol. As the study was retrospective and involved no more than minimal risks to the subjects, the need for patient consent was waived. Our study cohort included patients with pathologically-confirmed lung adenocarcinoma with locally advanced NSCLC (overall stage II-III). Patients treated with surgery or chemotherapy prior to the CT simulation date were excluded from our analysis. Patients receiving SBRT treatment were also excluded. For unbiased validation purposes, our 200 cohort was temporally split into two halves, the training Dataset A (n = 100) and the independent validation Dataset B (n = 100).

Clinical endpoints
The clinical outcome evaluated for this study was DM. Follow-up CT scans were performed every three to six months after treatment for tumor progression assessment. DM was considered as the disease spread to sites outside of the lungs. Time to DM was defined as the time interval between the start date of radiation therapy and the first scan date of radiographicallyevident DM and censored at the date of last negative scan in patients without recurrence. Time to OS was defined as the time between radiotherapy start date and date of death, and censored at the last follow-up date.

CT image acquisition and segmentation
Planning CT images were acquired on GE LightSpeed RT16 CT scanners (GE Medical Systems, Milwaukee, WI, USA) following clinical imaging protocol which mostly uses 120 kVp and reconstructed using standard convolution kernel. The most common voxel spacing on the CT image was 0.93mm, 0.93mm, 2.5mm. Exceptions to this protocol were two cases scanned using 140 kVp and 6 cases reconstructed using 3.75mm or 5mm slice thickness. All patients received contrast injection unless there existed contraindication. The primary tumors were contoured using Eclipse software (Varian Medical Systems, Palo Alto, CA, USA) by an experienced CT imaging researcher. Tumors contours were reviewed by an expert radiation oncologist (R.H.M). To ensure that the imaged tumor regions are of good quality for our analysis, cases with motion artifact were excluded.

Peritumoral contour preparation
Based on the proximity to the primary tumor, two peripheral regions of tissue were designed and termed here as tumor rim and tumor exterior. Tumor rim was defined to be the region that included the outer 3 mm of the tumor and 3 mm of tumor-adjacent lung tissue on either side of the tumor contour boundary; tumor exterior the region of lung tissues extending from 3mm to 9mm outside of the tumor contour (Fig 1.1). Rationale for tumor rim tissue inclusion was to designate a "real invasive front" and account both for the aggressive invading front of the tumor tissue and the adjacent normal lung layer, where cancerous islets can be frequently found [34]. The designation for tumor exterior followed from the spatial extent of microscopic tumor presence as found in pathology studies [28,36,37]. The generation of masks for delineating the tumor edge and the exterior tissue was accomplished through mathematical morphology operations of erosion and dilation, implemented using SimpleITK toolbox [38]. Since we were interested in the lung parenchymal tissue characteristics that may show association to cancer metastasis, care was taken to only include the tissue contours that are within the lung masks.

Radiomic feature extraction and selection
A comprehensive radiomics computational toolbox, PyRadiomics, was employed for feature extraction [39] (Fig 1.2). Designed to facilitate data reproducibility, PyRadiomics platform provided open-source standardized algorithm for radiomic feature computation. Its implementation included built-in wavelet and Laplacian of Gaussian filters for image processing and computed a total of 2175 radiomic features, including first order statistics, shape, and texture classes. Texture classes included gray level co-occurrence matrix (GLCM), gray level run length matrix (GLRLM), gray level size zone matrix (GLSZM), neighboring gray tone difference matrix (NGTDM), gray level dependence matrix (GLDM), and gray level distance zone matrix (GLDZM). Feature computation was performed at resampled voxel dimensions of 3×3×3mm 3 and an intensity bin width of 25 Hounsfield units. Radiomic feature extraction was performed for each of the three regions investigated in this study, i.e., tumor, tumor rim, and tumor exterior (data in S5-S10 Files). Feature selection followed a two-step procedure, feature stability and relevance, as shown in Fig 1.3. The selection of stable features was performed on Dataset A for each tissue region using the external test and retest RIDER dataset [40], subject to the intraclass correlation coefficient criterium (ICC > 0.85) [41] (description in S1 File). Subsequently, the set of stable features was processed using the minimal redundancy maximal relevance technique (mRMR) for dimension reduction, resulting in fifteen radiomic features for each region. Based on mutual information (MI), mRMR performed feature selection sequentially by determining the feature with maximum MI with the target variable and the minimum MI with the already selected features (Bioconductor "survcomp" package [42]).
The prognostic value of the peritumoral radiomic features was compared to the tumor ones as well as conventional and clinical parameters. The conventional variables considered for this study were the maximal 3D tumor diameter, and the volume of gross tumor volume (GTV). Clinical model included gender, age, overall stage, T-stage, N-stage, performance status, and tumor size.

Univariable analysis
The prognostic value of the radiomic features was evaluated using concordance index (CI) from "survcomp" package [43]. Noether's test was applied to assess the statistical significance of the computed CI from random chance (CI = 0.5) [42]. To account for multiple testing, a false-discovery-rate procedure by Benjamin and Hochberg was applied to adjust the p-values [44]. Univariable analysis was performed using the fifteen features selected using mRMR method for each of the tumor regions. All analyses were performed using R software (version 3.3.1) [45].

Multivariable model construction and validation
Multivariable models were constructed using Cox regression method, where the model was trained using Dataset A and the model predictions were validated in Dataset B. The multivariable radiomics models were constructed with the mRMR selected features where the feature complementarity was explored for potential prediction enhancement. Based on the principle of parsimony, the fifteen features for each tumor region were included to the model incrementally according to their mRMR ranking and 1000 cross-validations were performed on Dataset A for each intermediate model in order to assess its predictive power. The cross validations were performed through random subset sampling with balanced event ratios of 70:30 using caret package [46]. For each tumor region, the optimal feature set was the combination that rendered the highest mean CI value before decreasing and was termed signature (Fig A in S2 File). To determine the improvement due to radiomic features, clinical factors were incorporated to the individual radiomics model to generate combined clinical-radiomics model. The statistical significance of the model performance between a pair of multivariable models was assessed using the cindex.comp function ("survcomp" package). To assess clinical efficacy, we evaluated the statistical significance of our best performing combined clinical and peritumoral radiomic signature to multivariable clinical model and combined clinical and tumor radiomics signature.

Assessing patient stratification by peritumoral radiomics model
Finally, to demonstrate the clinical efficacy of our proposed peritumoral radiomics model, a Kaplan-Meier analysis was performed. We investigated the performance of three models, the clinical model, the combined clinical and tumor radiomics signature, and the combined clinical and rim radiomics signature, where the cox regression technique was used for combining separate models. For each model, the median prediction value from training Dataset A was used for stratifying the patients in Dataset B. A log rank test was performed to determine the statistical significance of risk to developing DM between the two groups.

Clinical characteristics
A total of 200 NSCLC patients with adenocarcinoma were analyzed in this study. At the time of diagnosis, the median age was 64 years (range: 35-93 years). The median follow-up was 28.1 months (range: 1.8-142.1 months). The median time to DM was 13.5 months (range: 0.3-119.1 months), with 128 (64%) patients having developed DM versus 72 (36%) who did not. Patient characteristics and cancer progression information can be found in Table 1. Peritumoral radiomics for metastasis prediction in locally advanced NSCLC Fig 2 displays the results of univariable analyses performed on Dataset A. In each tumor region, the full set of radiomic features was reduced to fifteen features based on relevance to DM events. In the tumor region (Fig 2A), ten features were significantly predictive of DM, had a CI range of 0.59-0.64, p-value < 0.05, and were all texture-based except for the 3D firstorder median, a statistics-type feature ( Table B in S3 File). The top tumor feature was GLCM Difference Entropy, which measures the variability in neighborhood intensity differences. In the rim region (Fig 2B), six significant features were found with CI range: 0.59-0.63 (Table D in S3  File). Four of these features were texture-based and two were statistics-based (first-order range and first-order minimum). The best performing feature in this region was GLRLM RunEntropy, providing a measure for the randomness in the distribution of run lengths and gray levels, with higher value indicating more heterogeneity in the texture pattern. In the exterior region, the only significantly predictive feature was first-order Kurtosis, which is a descriptor of the intensity distribution, with higher kurtosis value indicating the distribution is concentrated more towards the tails (Table F in S3 File). More detailed description of the selected features is provided (Tables A, C, E in S3 File). For comparison, the conventional features of tumor volume and maximum 2D axial and 3D diameters had CI values of 0.54, 0.55, and 0.54, respectively and were not statistically significant. In addition, the maximum 3D diameter and the contour volume of the rim and exterior regions were also not predictive. The feature expression trends between DM and non-DM cases were shown in Figs A-C in S4 File.

Radiomic signature validation and inter-comparison
Multivariable models were generated based on cox proportional hazard method. The forward selected radiomic signature from each tumor region was validated using Dataset B. The tumor radiomic signature was log sigma 0.5mm 3D GLCM DifferenceEntropy, which quantifies the intensity variability in neighboring voxels, and wavelet HLL GLRLM RunEntropy, which measures the variation in the distribution of run lengths and gray levels. The tumor rim radiomic signature consisted of LoG 1.5mm 3D GLRLM RunEntropy and Wavelet LHL NGTDM complexity; these measure, respectively, the entropy of gray level runs and large intensity changes in neighboring pixels. The radiomic signature of the exterior region includes log sigma 2.5mm 3D firstorder Kurtosis and log sigma 3.0.mm 3D NGTDM Strength, which, respectively, measures the fourth moment of the intensity distribution and deviation from homogeneity. To account for the potential confounding effect of tumor size and volume, statistical significance was tested between our radiomic signatures and these factors. The prognostic performance of the tumor and rim radiomic signatures were determined to be significantly stronger than tumor dimension or volume.
The performance of the radiomic signatures was compared to a clinical model constructed using Cox regression method. In Dataset B, this clinical model achieved a CI value of 0.53 (pvalue < 0.44). In the radiomics models, the multivariable rim signature achieved a CI value of 0.64 (p-value < 2.37×10 −5 ) in Dataset B, compared to the multivariable tumor signature CI value of 0.59 (p-value < 0.04) and the multivariable exterior model of CI value of 0.55 (pvalue < 0.15). Incorporating the rim multivariable model to the clinical parameters yielded a CI value of 0.65 (p-value < 7.57×10 −6 ). For comparison, this combined model was found to be significantly more predictive than the clinical model (p-value < 0.003). A composite radiomics model including the tumor, rim and exterior regions was also constructed (CI = 0.63), but was found to be less predictive than the tumor rim signature (p-value < 0.30). The prediction by combined clinical and tumor radiomics signature was found to be not statistically different from the combined clinical and rim predictor (p-value < 0.13), neither was the composite clinical, tumor and rim predictor (p-value < 0.38). The results of multivariable model validation and inter-model comparison were displayed in    (Fig 4A). Neither the log rank test showed a significant p-value nor did the hazard ratio show significance in the validation cohort. Fig 4B showed risk stratification using the combined clinical and tumor radiomics model, which had a significant log rank p-value of 0.012. Fig 4C  demonstrated the potential stratification that can be achieved using our proposed model combining clinical and rim radiomic signature. Using the patient risk scores derived from the training cohort, Kaplan-Meier analysis in validation Dataset B showed statistically significant difference (p-value < 1×10 −3 ) for metastasis-free probability estimates. The lower risk group showed a hazard ratio of 0.44 compared to the higher risk group. Comparison of prognostic performance across the different multivariable models in the validation cohort (n = 100). � indicates statistical significance (" � " indicates p-value <0.05, " ��� " indicates p-value <0.0001 from random prediction (Noether test)). From left to right, the compared multivariable models include clinical, visible tumor, exterior, tumor rim, combined clinical and tumor, combined clinical and rim, and the combined radiomics model. Crossbars indicate the comparison made between CI of two multivariable models, where � indicates significant difference and ns not statistically significant. https://doi.org/10.1371/journal.pone.0206108.g003

Discussion
Advances in cancer biology research have provided insights on how the tumor proliferates through its interaction with the surrounding normal tissue [47][48][49][50]. Tumor invasion into the peripheral normal tissues on the cellular level may translate to tissue morphological changes, which, in turn, may inform us about the level of metastatic activity. Here, we leveraged radiomics analysis to quantitatively characterize the peritumoral imaging features on planning CT images and assessed their prognostic power on DM. To our knowledge, this is the first study that correlates radiomic features from the normal-appearing peritumoral tissues with cancer metastasis in NSCLC.
The use of radiomic features to predict DM for NSCLC has only been investigated in few studies [24,25]. Fried et al [25] extracted 198 tumor features from averaged CT, 4DCT, and contrast-enhanced CT images; and showed qualitatively through Kaplan-Meier analysis that the combination of radiomic features with clinical parameters may enable patient stratification for DM in a 91-patient stage III cohort. Coroller et al [24] investigated intratumoral radiomic features for DM in a cohort of 182 NSCLC and demonstrated predictive power in a validation cohort (CI = 0.61). However, these were based on tumor-only radiomic features and did not account for the cancerous infiltration into the surrounding normal parenchymal tissues. This study explored the tumor rim and exterior radiomic features as potential DM prognosticator in a larger cohort of 200 patients with locally advanced NSCLC using a quantitative metric of CI. Importantly, we discovered a higher prognostic value for the tumor rim radiomic signature than the tumor-only one, where a comparison between the two showed statistical significance (p-value = 0.048).
Our data driven approach led us to the discovery of potentially important tumor rim radiomics signature for DM prediction and also the finding that the other tumor regions show relatively less prognostic power. In the exterior region, that the multivariable radiomics model was not significant was consistent with the overall weaker signal as evident in the univariable analysis. As for the visible tumor region, its significant features were similar to the rim region in terms of feature type and prognostic power. However, while the selected features from the tumor showed similar performance as the rim features univariablely, the rim radiomic signature showed the superior performance over the tumor one. This may be due to the complementary effect accomplished by combining top performing feature, LoG 1.5mm 3D GLRLM RunEntropy, and the moderately performing feature, Wavelet LHL NGTDM complexity. Interestingly, as larger value from either rim feature suggests reduced risk to developing DM, i.e. tumor rims with more heterogeneity in run lengths and gray levels as well as more rapid changes in gray level intensity would be less prone to tumor metastatic activities. Moreover, despite the fact that the CI of the combined clinical and rim radiomics model prediction was not shown to be statistically different from that of the combined clinical and tumor radiomics one, Kaplan-Meier analysis suggested that the former would allow for better risk stratification in terms of log rank p-value (Fig 4B vs. 4C). Furthermore, our identification of the tumor rim as the crucial imaging biomarker for DM was consistent with pathology and tumor biology findings and it was well known that the periphery of the tumor harbored many activities of cancer invasion and metastasis, e.g. epithelial-mesenchymal transition [51], tumor-associated macrophages [50,52], tumor budding [53], and lymphovascular invasion [31,54]. While a correlation with biological processes at the tumor periphery was beyond the scope of the present study, ample evidence from tumor biology literature supported our hypothesis that the tissue features on the tumor-normal interface may indicate tumor aggressiveness towards DM. Thus, our findings were hypothesis generating and may facilitate new discoveries in the DM prognostication using information originated from tumor rim.
Limitations of this study include the choice of 6mm shell of tissues around the tumor for our correlation analysis with DM. Given that there may exist individual variation in terms of the disease spread pattern and location of extratumoral cancer colonies, we attempted to address this by taking a relatively wide margin of 6mm and expand our region of interest radially outward from the tumor to mimic the disease spread pattern. Other limiting factor may be the variation in CT acquisition parameters as patient CT simulation dates spanned from 2001 to 2014. We sought to mitigate this by removing cases of motion artifacts and performed image resampling at 3×3×3 mm 3 to reduce voxel noise. Lastly, our findings may be limited by our cohort size (n = 200) and patient cases collected in our institution. We had performed temporal split of our data to generate an independent validation cohort, of similar patient and treatment characteristics, for model testing. Future investigations would involve testing our hypothesis by expanding our study to include patients of other histology types and evaluating the generalizability of our findings using multi-institutional image data. In spite of these limitations, our investigation demonstrated differential predictiveness of imaging features between the tumor and its surrounding tissues for distant metastatic spread.

Conclusion
In conclusion, we have demonstrated strong prognostic value of peritumoral radiomic features for DM in patients with locally advanced NSCLC. The presented rim radiomic signature was independently validated and was shown to have better predictive power compared to tumor radiomic signature. Such pretreatment imaging predictor may benefit patients susceptible to developing DM in precision medicine approach.
Supporting information S1 File. Fig A. Tumor regions generated for RIDER test/retest data. Subfigures a., c., and e. represent the contours of the tumor, tumor rim, and the tumor exterior, respectively, from the test dataset. Subfigures, b., d., and f. represent the counterparts in in the re-test dataset.