F-fluorodeoxyglucose positron-emission tomography (FDG-PET)-Radiomics of metastatic lymph nodes and primary tumor in non-small cell lung cancer (NSCLC) – A prospective externally validated study

1 Department of Radiation Oncology (MAASTRO), GROW–School for Oncology and Developmental Biology, Maastricht University Medical Center (MUMC +), Maastricht, the Netherlands, 2 Institute of Radiooncology—OncoRay, Helmholtz-Zentrum Dresden-Rossendorf, Dresden, Germany, 3 Department of Radiotherapy and Radiation Oncology, Medical Faculty and University Hospital Carl Gustav Carus of Technische Universität Dresden, Dresden, Germany, 4 OncoRay, National Centre for Radiation Research in Oncology, Medical Faculty and University Hospital Carl Gustav Carus of Technische Universität Dresden, Helmholtz-Zentrum Dresden-Rossendorf, Dresden, Germany, 5 Department of Radiology and Nuclear Medicine, Radboud UMC, Nijmegen, the Netherlands, 6 Department of Radiology, Leiden University Medical Center, Leiden, the Netherlands, 7 Biomedical Photonic Imaging Group, MIRA Institute, University of Twente, Enschede, the Netherlands, 8 Department of Radiation Oncology, Radboud University Medical Center, Nijmegen, the Netherlands

Introduction Non-small cell lung cancer (NSCLC) patients often present with hilar and/or mediastinal lymph node involvement at diagnosis or during the course of disease. Lymph node stage prior to treatment is strongly related to disease progression and worse prognosis [1]. Furthermore, it affects treatment selection and target volume definition, for metastatic lymph nodes in patients eligible for high-dose (chemo)radiotherapy [2].
In this study we hypothesized that the local selection of more aggressive cancer cells in the metastatic hilar/mediastinal lymph nodes, being likely to determine prognosis, may provide an additional and valuable source of information to the primary tumor for NSCLC patients. A Radiomics-based approach comprises the extraction of a large set of imaging descriptors [3]. The underlying hypothesis is that biomarkers of imaging phenotypes deliver complementary and clinically relevant information, which could be incorporated into individualized radiation oncology approaches and shared decision-making tools [4][5][6][7]. To demonstrate this, we performed a combined Positron Emission Tomography (PET) Radiomics analysis of metabolic activity as measured with 18 F-fluorodeoxyglucose (FDG) uptake in both primary tumor and metastatic lymph nodes, and further validated these results in an independent cohort.

Development cohort
Patient population. The prospective data collection was approved by the Institutional Review Board of the Department of Radiotherapy of Maastricht University Medical Center (MAASTRO clinic) (clinicaltrials.gov NCT00522639). Electronic medical charts of NSCLC patients were reviewed. Patients undergoing surgery, Stereotactic Body Radiotherapy (SBRT) or palliative treatment, or who had a previous malignancy within five years prior to diagnosis were excluded from analysis. A total of 343 NSCLC patients (stage I-IIIB) referred to curative treatment (between May 2006 and September 2012) were selected for the development cohort. Out of these, 262 patients (76%) had metastatic lymph nodes. Patients received high-dose radiotherapy (RT), planned on a dedicated FDG-PET-CT scan, combined with chemotherapy. Clinical follow-up was performed according to national guidelines. All patients in the development cohort provided informed written or verbal consent to data inclusion in clinical studies.
Image acquisition. Before scanning, patients fasted for at least 6 hours. Two different protocols were used: until December 2010, the total dose of FDG was calculated as (bodyweight x 4 + 20) MBq, and as from January 2011, the administered dose was (2.5 x bodyweight) MBq as defined by the NEDPAS protocol [8]. FDG-PET-CT images were taken 60-minutes post injection. Data acquired until December 2006 were gathered on Siemens Biograph 16 CT-PET scanner, and from that time onwards on a Siemens Truepoint 40 CT-PET (Siemens Healthcare AG, Erlangen, Germany). An Ordered Subset Expectation Maximization 2D 4 iterations 8 subsets (OSEM2D 4i8s) algorithm was used for image reconstruction using post-reconstruction 5mm Gaussian filtering, and voxel size of 4.0728 x 4.0728 x 3 (mm). Model-based methods were applied for scatter correction. All PET scans were corrected for attenuation using the mid-ventilation phase of the 4DCT or a 3DCT thorax in case the 4DCT was not of sufficient image quality due to irregular breathing of the patient. All exams were corrected for random events and decay.

Validation cohort
Patient population. The validation cohort included 215 stage I-IIIB NSCLC patients, treated with primary radio(chemo)therapy between May 2006 and October 2012 at the Department of Radiotherapy of Radboud UMC Nijmegen, following same inclusion and exclusion criteria as for development cohort. In total 115 (53%) patients were node-positive, of which 50 (23%) had an available treatment planning FDG-PET-CT. All patients in the validation cohort provided informed written or verbal consent to data inclusion in clinical studies.
Image acquisition. Before scanning, patients fasted for at least 6 hours. FDG-PET scans were performed 60 minutes after intravenous injection of approximately 250 MBq FDG (Covidien, Petten, the Netherlands) and 10 mg furosemide. PET scans were performed on Siemens Biograph Duo (Siemens Medical Solutions USA, Inc.) using three-dimensional emissions of 4 minutes per bed position as described previously [9]. A low-dose CT scan for localization and attenuation-correction purposes was acquired. Scanning parameters included 40 mAÁs (50 mAÁs for patient weight >100 kg and 60 mAÁs if >120 kg), 130 kV, 5-mm slice collimation, 0.8-second rotation time, and pitch of 1.5, reconstructed with 3-mm slices for smooth coronal representation. An OSEM2D 4i16s algorithm was used for PET image reconstruction, with a voxel size of 5.3 x 5.3 x 3.375 (mm). All PET scans were corrected for attenuation using CT and simulation approaches. Model-based methods were applied for scatter correction. All exams were corrected for random events and decay.

Regions of interest (ROI)
Images were imported into the research treatment planning system Xio/Focal (development cohort) and Eclipse (validation cohort) using the DICOM protocol. The primary gross tumor volume (tumor) and metastatic hilar/mediastinal lymph nodes (LN), identified as PET positive and/or proven by endoscopic ultrasound bronchoscopy/esophagoscopy (EBUS/EUS), were manually delineated by experienced radiation oncologists on the fused FDG-PET-CT images, and used as the regions of interest for analysis [10][11][12]. A single structure representing all metastatic lymph nodes, regardless of their number, was derived for each patient (see online appendix for further details).

Image analysis
In-house developed software was used to extract Radiomics descriptors from the FDG-PET scans [13][14][15][16]. Imaging descriptors comprised first order statistics (n = 16), shape and size (n = 13), intensity volume histograms (n = 45), and textural features describing the spatial distribution of voxel intensities (n = 44). Textural features were calculated from grey-level cooccurrence (GLCM), grey-level run-length (GLRLM) and grey-level size-zone texture matrices (GLSZM). To determine these matrices, images were first discretized with a bin width of 0.5 (Standardized Uptake Value or SUV), according to: Where I is the original image, I(x) represents the SUV of voxel x, and I D is the resulting discretized image [15]. Texture matrices were then constructed by considering 26 connected voxels (i.e. voxels were considered to be neighbors in all 13 directions in three dimensions) at a distance of 1 voxel. Features derived from GLCM and GLRLM were calculated by averaging their value over all 13 directions. Forty-four textural features were extracted (22 GLCM, 11 GLRLM and 11 GLSZM). In total, 118 imaging features were calculated based on the FDG-PET distribution within ROI, which mathematical formulations are detailed in the work of Leijenaar et al. [14]. Image analysis was performed in Matlab R2012b (The Mathworks, Natick, MA), based on an adapted version of Computational Environment for Radiotherapy Research (CERR) [17].

Statistical analysis
Study parameters/endpoints. Primary endpoint of the study was overall survival (OS), defined as the time from the start of radiotherapy until the last day of follow-up or death due to any cause, and was available for all patients under analysis. A patient still alive at the end of the study was regarded as right-censored.
Univariable analysis. A pre-feature selection was performed as detailed in the online appendix. Imaging and clinical features were analyzed as continuous variables in a univariable Cox regression. In addition, a correlation analysis was conducted for maximum, peak, and mean SUV, and volume of tumor and LNs, and tumor load (sum of tumor and LN volumes).
Multivariable analysis. A prognostic model was fitted to the data with a 10-fold cross-validated least absolute shrinkage and selection operator (LASSO), selecting amongst PET-Radiomics descriptors extracted from the tumor (model 1), LN (model 2) and the union of both structures (model 3) [18]. LASSO selects variables correlated to the measured outcome by shrinking down to zero coefficients weights for features non-related to outcome. Features were entered in the model as continuous variables. Regression coefficients, hazard ratios (HR) and confidence intervals (CI) were estimated using the whole development cohort. A diagram illustrating this methodology is shown in Fig 1. Log-likelihood tests for non-nested models, Akaike information criterion (AIC), were performed to compare the fit of the three derived and independent models. AIC measures the relative quality of model fit to a given cohort, providing substantiation for model selection, i.e. the one with the lower AIC is the preferred model [19]. Log-linearity assumption was verified for the selected features in the final models by fitting a penalized smoothing spline. Cox proportional hazards assumption was graphically examined with the Schoenfeld residuals.
Model performance. Model performance was assessed in the development and validation cohorts by means of a concordance-index and corresponding 95% CI [20]. Concordance-index or Harrell's C-index, evaluates the fraction of patient pairs for which the predicted and actual outcome are concordant, ranging from a random 0.5 to a perfect 1 [21].
All statistical analysis was conducted in R (version 2.15.2), using the libraries: survival, survcomp, glmnet, cvTols and rms [22]. The development cohort is publicly available at at www. cancerdata.org and Dryad (doi:10.5061/dryad.752153b). Table 1 gives a complete overview of patients under analysis including treatment details. Node-positive patients with available FDG-PET-CT scans were included in the analysis: 262 for development and 50 for validation phase. A univariable analysis was performed for each of the clinical variables in the development cohort (Table 2). Results show that TNM staging was not correlated to OS, while N stage was a prognostic factor at 0.1 level (p = 0.09) in the development cohort, with higher stages being associated with a worse prognosis (hazard ratio (HR) stage 2 = 1.44; HR stage 3 = 1.75). The number of LN stations was significantly associated with a higher risk. Radiotherapy dose was also significantly associated with prognosis information. None of the remaining analyzed metrics, including gender, age, histology and chemotherapy showed a significant correlation to OS in our development cohort.

Results
Results of univariable Cox regression of FDG-PET Radiomics features extracted from both tumor and LNs are shown in the appendix (S2 File. Tables). Only short run emphasis from the texture GLRLM group was significantly correlated to OS when extracted from the primary tumor in the development set. On the other hand, metrics derived from LN showed a good univariable correlation to outcome, with 13 of the 16 analyzed features being significantly related to overall survival. Table 3 displays partial results of this analysis, for which a high Pearson correlation between metabolic features within each structure could be verified, but neither a correlation could be found with own volume, nor with metabolic features of the other structure (Fig 2).
Three model approaches were derived and fitted to the data in the development cohort. These are represented in Table 4, with corresponding hazard ratios and concordance-index. Selected features were log-linear, except for LN volume, which had to be converted into a logarithmic scale. The proportional hazards assumption was satisfied for all features. Graphical assessment of these assumptions can be appreciated from

Discussion
Disease management of NSCLC is a primary concern, for which prognostic assessment is essential to fulfil the potential of individualized and personalized treatment. Nowadays, a wide range of information sources are available of which the non-invasive type play a fundamental role in reducing the patients' burden [4][5][6]. Among these, metabolic imaging has been increasingly explored for prognosis assessment, based on SUV patterns of FDG of the primary tumor, as an extension of its primary diagnostic function: detection of metastatic lymph nodes and distant metastasis [23,24].
NSCLC patients often present with lymph node involvement at diagnosis, which deeply impacts prognosis and response to treatment [1,2]. Apart from the number of metastatic lymph node stations, lymph node size and corresponding metabolic activity may vary among patients [25]. Diagnosis of metastatic hilar/mediastinal lymph nodes is commonly performed through FDG-PET-CT and consecutive EBUS/EUS or mediastinoscopy [26][27][28][29][30]. Diffusionweighted magnetic resonance imaging (DWI) is also a solid alternative for lymph node detection, as proven by the work of Shen et. al [31]. In the work carried out, given the availability of FDG-PET-CT scans for all patients under analysis as part of their protocol for radiotherapy treatment, we proceeded the analysis with lymph nodes detected using this imaging modality. Based on the rationale that disease progression and the ability to metastasize are closely related to the presence of metastatic lymph nodes, we hypothesized that FDG-PET-based Radiomics information of these nodes would provide additional prognostic information to the information that is obtained from the primary tumor [23]. Radiomics has been proven to have prognostic potential in predicting clinical outcomes or treatment monitoring in different cancer types [13,16,32,33] and can essentially be applied to different medical imaging modalities and disease-related structures such as the primary tumor, metastatic lymph nodes or metastatic lesions.
Common clinical metrics associated with disease prognosis, including TNM staging (p = 0.92), could not be attributed statistical significance in the development cohort. Nevertheless, patients were staged in accordance to later editions of the classification system, that newer classification as updated by the 8 th edition may overcome [34]. On the other hand, higher number of metastatic lymph nodes could be associated with a statistical significant higher risk (p<0.01), already anticipating the extra information these structures may provide, that we further evaluated through imaging descriptors. Given the broad range of imaging descriptors analyzed by a Radiomics approach, an initial step in this analysis was to perform an exploratory univariable analysis of the most commonly analyzed PET metrics when extracted from tumor and LNs, namely maximum, peak, and mean SUV, volume and tumor load. None of the metabolic metrics extracted from tumor had significant prognostic value, whereas the same ones extracted from LNs were related to OS and yielded an univariable C-index of at least 0.57 (Table 3). No statistical significance could be associated with tumor volume (p = 0.47), as opposed to LN volume (p<0.01; c-index 0.60, 95% CI = 0.51-0.68). Tumor load, sum of tumor and LNs volumes, also had prognostic value, in line with previous studies [35]. Despite the strong correlation between metabolic features within each structure, neither a correlation with own volume, nor with metabolic features of the other structure could be found, contrary to previous evidence of correlation between maximum SUV of tumor and total LN volume [36]. One can already infer from this univariable analysis the considerable prognostic value that the metabolic distribution within metastatic lymph nodes has for node-positive NSCLC patients, which is further reinforced by remaining features.
The main purpose of this study was to analyze imaging features extracted from the LNs and compare it against the ones extracted from the primary tumor. In addition the non-univariable statistical significance attributed to common clinical features (e.g. TNM staging), subsequent modeling conducted solely on FDG-PET imaging features. Prior to development of the three prognostic models, a pre-selection of the imaging features, comprising their stability and robustness as discussed by Leijenaar et al. [14], and supported by Desseroit et al. [37], was conducted and is discussed in the appendix section. Of the three derived independent fits, the model with features from LNs alone or LNs in combination with primary tumor performed best (c-index 0.62). The model included at least one feature from categories of shape and size, first order statistics and texture descriptor, excluding features from the IVH group, despite Analysis was conducted for primary tumor and metastatic lymph nodes separately, and for both structures in combination. Hazard Ratios (HR) and corresponding pvalues are reported. Performance of the model is expressed by internal and external ÃÃ concordance-index (C-index). Internal performance includes associated 95% confidence-interval (CI) of the C-index.
Acronyms: GLCM-Grey Level Co-occurrence matrices; GLRLM-Grey Level Run-length matrices; Stats-first order statistics Ã A logarithmic transformation was applied to LN volume ÃÃ External validation https://doi.org/10.1371/journal.pone.0192859.t004 FDG-PET-radiomics of tumor and lymph nodes in NSCLC their proven univariable prognostic evidence in a previous study [38]. The intrinsic heterogeneity described by the FDG distribution within a primary tumor has proven its prognostic power, not only for NSCLC, but also for other cancer types [39,40]. Inclusion of such descriptors in derived prognostic models, gathering heterogeneity insight into not only primary tumor but also metastatic lymph nodes, resulted in a model with a better fit and more accurate description of disease structures. Shape descriptors, including total LN volume, included in the multivariable model, were shown to be an independent prognostic parameter. In a previous study, tumor load was revealed to be a prognostic factor, but not LN volume independently, which we proved in this study [35]. Finally, and despite their univariate prognostic value, particularly when assessed from LN, none of the most common SUV descriptors were included in the final models. However, as Radiomics analysis includes a large number of features, selection of the most promising ones is difficult. We attempted to overcome this difficulty using an automatic feature selection routine (LASSO), to ensure an optimal exclusion of redundant or highly correlated features from the final models derived, further complemented with and external validation. Nonetheless, it cannot be completely excluded that other variables have similar or even greater prognostic value than the current ones, and therefore larger imaging datasets are needed to validate and confirm our findings. Similarly, another limitation of our study was the lack of additional patient data for the validation phase. A lower performance is commonly observed when validating a model against new, independent, and external cohorts, which is most frequently attributed to discrepancies between development and validation data [41]. A larger validation cohort would increase the robustness of the validation procedure. Nevertheless, we observed a benefit from combining imaging features from both primary tumor and metastatic lymph nodes for node-positive NSCLC patients. In terms of prognostics, this should not be disregarded, particularly when compared to the limited capacity of humans to infer an accurate prognosis from same data [42]. In fact, the most commonly accepted standard for disease prognosis assessment failed to provide a significant stratification of risk patients (TNM stage). A recent study with 139 NSCLC patients treated with at least 60Gy with a concurrent (chemo)radiotherapy regimen showed the importance of measuring the post-treatment SUV in the metastatic lymph nodes, as an increase in both the absolute value and percentage of residual activity compared to the pre-treatment SUV were associated with worse local-regional control [43]. In our study we looked only at the pre-treatment scan and could already discern an association with OS of imaging features derived from the metastatic lymph nodes that could not be assessed based on the primary tumor for node-positive patients. However, it is hypothesized that the the variation of FDG-PET-Radiomics features between subsequent scans at an early phase of treatment for both the primary tumor and the metastatic lymph nodes, and their impact on survival for NSCLC patients as a complement to the positive findings reported here. Also, a positive correlation between PET information derived from LNs and overall relapse has been reported [44]. Our analysis focused primarily on overall survival, and therefore we could not validate these findings. In the future we will be able to analyze this outcome as we are currently improving our data collection routines, to further evaluate other outcomes, actually limited to overall survival. Nevertheless, our findings emphasize the importance of analyzing FDG-PET signal of metastatic lymph nodes prior to radiotherapy, to further complement the information retrieved from the primary tumor.
Based on the positive and relevant findings we documented, we have plans to extend our analysis in a similar manner to other disease sites, particularly head and neck cancer, for which the involvement of the nodes is a well-known prognostic risk factor [45]. Likewise, we are also aiming to analyze other PET tracers, particularly hypoxia markers (e.g. HX4) and compare the results with FDG uptake [46,47]. Further, it has been recently shown that combining Radiomics from both FDG-PET and low-dose CT can improve prognosis in NSCLC [48], therefore, combining Radiomics from LNs extracted from the CT part of the FDG-PET-CT would be a much desirable and natural step for future analysis with larger datasets. Seemingly, other endpoints could obviously be analyzed following this Radiomics approach, namely FDG-PET-Radiomics signal correlation to tumor marker. This rationale follows the correlation demonstrated between maximum SUV and an increased expression of glucose transporters 1 and p53 for adenocarcinoma NSCLC patients, but not squamous cell by the work of Taylor et al. [49]. As pointed out, these studies need to take into consideration the effect of different tumor histology types, which correction was not performed in this work. Resistance of p53-related chemotherapy has also been shown to be linked with maximum SUV by the work of Duan et al. [50]. Correlation to therapy response was however not the subject of this study, as the endpoint under analysis was overall survival due to any cause, but one can definitely leave for future approaches the link between Radiomics descriptors and tumor markers related to response to therapy.
In summary, common SUV descriptors derived from metastatic lymph nodes were associated with overall survival in a large cohort of NSCLC patients. Additionally, PET information demonstrated to have higher prognostic value when extracted from metastatic lymph nodes in comparison to the primary tumor alone, further complementing its information. The use of 3D information based on imaging is becoming a broader field with expected great gain for patients' outcomes assessment and treatment planning adaptations, following its application to dissimilar structures as primary tumor, metastatic lymph nodes and possibly distant metastasis.
Supporting information S1 File. Feature pre selection. Methodology conducted for feature pre-selection for both primary tumor (tumor) and metastatic lymph nodes (LN), and corresponding results and interpretation. (DOCX) S2 File. Tables. Univariable analysis of stable and robust Radiomics features from primary tumor for development and validation datasets and that were common surrogates as extracted from the merged and largest node as well as the merged and most active node, and therefore entered as continuous variables in the multivariable model building. (DOCX)