Radiogenomic analysis of hypoxia pathway is predictive of overall survival in Glioblastoma

Hypoxia, a characteristic trait of Glioblastoma (GBM), is known to cause resistance to chemo-radiation treatment and is linked with poor survival. There is hence an urgent need to non-invasively characterize tumor hypoxia to improve GBM management. We hypothesized that (a) radiomic texture descriptors can capture tumor heterogeneity manifested as a result of molecular variations in tumor hypoxia, on routine treatment naïve MRI, and (b) these imaging based texture surrogate markers of hypoxia can discriminate GBM patients as short-term (STS), mid-term (MTS), and long-term survivors (LTS). 115 studies (33 STS, 41 MTS, 41 LTS) with gadolinium-enhanced T1-weighted MRI (Gd-T1w) and T2-weighted (T2w) and FLAIR MRI protocols and the corresponding RNA sequences were obtained. After expert segmentation of necrotic, enhancing, and edematous/nonenhancing tumor regions for every study, 30 radiomic texture descriptors were extracted from every region across every MRI protocol. Using the expression profile of 21 hypoxia-associated genes, a hypoxia enrichment score (HES) was obtained for the training cohort of 85 cases. Mutual information score was used to identify a subset of radiomic features that were most informative of HES within 3-fold cross-validation to categorize studies as STS, MTS, and LTS. When validated on an additional cohort of 30 studies (11 STS, 9 MTS, 10 LTS), our results revealed that the most discriminative features of HES were also able to distinguish STS from LTS (p = 0.003).

SCientifiC RePORTS | (2018) 8:7 | DOI: 10.1038/s41598-017-18310-0 hypoxia, treatment resistance, and metastasis of GBM 7 . The ability to capture tumor hypoxia will have significant clinical implications towards personalized treatment options for GBM patients and improve our understanding of tumor behavior as well as patient's outcome and response to specific treatments 8 .
As a standard-of-care protocol for brain tumor characterization, magnetic resonance imaging (MRI) is capable of capturing a diverse spectrum of tumor phenotypes. For instance, enhancement on Gd-T1w MRI is known to be correlated with blood brain barrier (BBB) disruption, while T2w/FLAIR abnormalities are known to capture proliferative tumor margins and vasogenic edema 9 . This suggests that the phenotypic differences at the cellular level are perhaps also reflected on MRI. Thus even though the visual appearance of different tumor phenotypes on MRI are similar, there nonetheless might be subtle sub-visual cues reflective of the differences in the micro-architectural appearance embedded in routine MRI that might enable distinction of molecularly distict GBM phenotypes. In this work, we hypothesized that phenotypic changes due to hypoxia-specific events (as observed on mRNA data), are also manifested on routinely acquired MRI (Gd-T1w, T2w, FLAIR); and can be captured using radiomic (high throughput computer extracted) features.
Recently, several studies have begun to explore the role of radiomic features on routine MRI scans (Gd-T1w, T2w, FLAIR) in capturing the underlying tumor pathology and molecular heterogeneity 10,11 . Radiomic features allow capture of quantitative imaging measurements by computing local macro-and micro-scale morphological changes in texture patterns (e.g. roughness, image homogeneity, regularity and edges) within the lesion. Many of these features, such as gray level co-occurrence matrix (GLCM)-based features, quantify enhancement heterogeneity, which has been shown to predict aggressive growth, unfavorable prognosis, and poor treatment response 12 . These texture-based radiomic MRI phenotypes have been shown to serve as surrogate markers for characterizing different molecular aberrations (including mutational status of IDH, 1p19q, MGMT) towards understanding GBM behavior and patient prognosis (also known as radiogenomics) [13][14][15] . However, to our knowledge, none of the existing studies have explored an explicit radiogenomic link between radiomic features obtained from routine MRI scans (i.e. Gd-T1w, T2w, FLAIR) and the hypoxia pathway and their role in predicting patient prognosis in GBM tumors. We present a radiogenomic approach to identify radiomic surrogate markers specific to the hypoxia pathway obtained from routinely acquired MRI scans (Gd-T1w, T2w, and FLAIR). We specifically focus on radiomic interrogation to capture the hypoxic events that lead to intratumoral molecular heterogeneity as manifested within the tumoral (necrosis, enhancing tumor) and peritumoral regions (non enhancing tumor and edema). Our work is motivated by previous studies that have demonstrated association of VEGF, a key mediator of tumor angiogenesis in the hypoxia pathway, with edema and tumor burden 16,17 .
In this work, we have two objectives. Firstly, we will identify a set of radiomic features obtained from Gd-T1w, T2w, FLAIR MRI scans, that were most discriminative of the extent of hypoxia in the tumor microenvironment, as measured using a hypoxia enrichment score (HES) using Single-sample Gene Set Enrichment Analysis (ssG-SEA) 18 . Secondly, given that tumor hypoxia is known to have implications in GBM survival 6 , we seek to investigate the role of radiomic surrogate markers of hypoxia (as reflected on the HES), in discriminating patients' overall survival (OS). The patients will be categorized based on their OS, as short-term survivors (OS < 7-months), mid-term survivors (7 months < OS < 16 months), and long-term survivors (OS > 16-months). Recent studies, including our own work, have independently investigated the use of radiomic descriptors in distinguishing short-term versus long-term survivors 12 . However, unlike previous studies, we present the first attempt at performing tumor and peritumoral interrogation towards (a) identifying radiomic surrogate markers of tumor hypoxia on pre-treatment MRI scans, and (b) evaluating their role in discriminating short-term, mid-term and long-term survivors of GBM. Our approach is intended to form a precursor to building novel image-based prognostic as well as predictive surrogate markers for personalizing treatment management in GBM, by reliably stratifying patients based on their hypoxia profile and overall survival.
The rest of the paper is organized as follows. Section 2 discusses the previous work and novel contributions. In Section 3, we provide methodological details of this work. Experimental results are presented in Section 4. We discuss the results in Section 5 and provide concluding remarks in Section 6.

Previous work and Overview
Recently, there has been some work in identifying radio-genomic associations of hypoxia in GBM as well as other tumors. For instance, Yopp et al 19 . have demonstrated inverse correlation of dynamic contrast-enhanced MR features with severity of hypoxia in hepatic cancers. Similarly, in a study by Diehn et al., proliferation and hypoxia gene expression patterns were found to be associated with the volume of mass effect and tumor contrast enhancement on different MRI protocols (Gd-T1w, T2w, FLAIR) 20 . However, to our knowledge, none of the works in the existing literature have attempted to establish a relationship between radiomic MR features obtained from different tumor sub-compartments and tumor hypoxia and then further employed these radiomic MRI surrogate markers of hypoxia as potential surrogate markers of overall survival in GBM. Figure 1 illustrates an overview of our framework. In Module 1, different MRI protocols are aligned in the same frame of reference; T2w, and FLAIR MRI were registered to Gd-T1w MRI in our case. In Module 2, the tumor sub-compartments including necrosis, edema, and enhancing tumor, are manually segmented by an expert using Gd-T1w, T2w, and FLAIR sequences, for each study. Module 3 involves computing radiomic features such as directional gradients (Gabor) and local intensity statistics (Haralick, Laws) from each of the tumor-specific sub-compartments (necrotic, enhancing tumor and peritumoral [edema and non enhancing tumor] regions) across the 3 MRI sequences. Within Module 3, top compartment-specific radiomic features are selected using mutual information 21 , based on their association with the hypoxia enrichment score, as obtained from the expression profile of 21 genes implicated in the hypoxia pathway (obtained via corresponding mRNA expression data) 20 . The top radiomic features that have the highest mutual association with hypoxia are then further employed to distinguish patients with short-term, mid-term, and long-term survival, using Kaplan-Meier survival analysis.

Methods
Study Population. Our cohort consisted of a total of 180 retrospectively analyzed treatment-naÃ¯ve multi-parametric MRI scans from the Cancer Imaging Archive (TCIA) 22 . TCIA is an open archive of cancer-specific medical images and associated clinical metadata which was curated by the National Cancer Institute (NCI) in association with participating institutions from the United States. Our inclusion criteria consisted of the following: (1) availability of all 3 routine MRI sequences (Gd-T1w, T2w, FLAIR) as well as RNA sequences for treatment-naïve patients, (2) MRI scans with diagnostic image quality (excluding the studies with image artifacts, as assessed by the expert reader), and (3) availability of individual overall survival information. A total of 65 cases were excluded either due to the absence of baseline scans, MRI artifacts, or unavailability of corresponding RNA sequence data or when one of the 3 MRI protocols (Gd-T1w, T2w, FLAIR) were not available. A total of 115 studies were used for further analysis. The corresponding RNA Seq for each of the 115 patients were collected from Broad Institute 23 . The 115 GBM subjects (69 males, age: 58.4 ± 12.56 yrs and 46 females, age: 56.47 ± 16.76 yrs), were then divided into three groups where the overall survival of the patients was stratified into short-term (≤7 months), long-term (>16 months) 24 , and the remaining into medium-term (>7 months to 16 months) survivors. Table 1 shows the patient demographics including mean overall survival, age, Karnofsky Performance Score (KPS) and gender of the population used in the study.  Preprocessing. For every study, T2w and FLAIR were co-registered with reference to Gd-T1w MRI using 3D affine registration with 12 degrees of freedom encoding rotation, translation, sheer, and scale. The registration was performed using the General Registration (BRIANSFIT) module of 3D Slicer 4.5 25,26 . To resolve the issue of resolution variability, every MRI slice within a scan was re-sampled to have a uniform pixel spacing of 0.5 × 0.5 mm 2 and was then interpolated to have 3 mm slice thickness. Skull stripping was done using the skull-stripping module in 3D Slicer 27 . Additional details regarding registration parameters employed in this work are provided in the Supplementary document. We then corrected the MRI protocols for known acquisition based intensity artifacts; bias field inhomogeneity and intensity non-standardness. Intensity non-standardness refers to the issue of MR image "intensity drift" across different imaging acquisitions. Intensity non-standardness results in MR image intensities lacking tissue-specific numeric meaning within the same MRI protocol, for the same body region, or for images of the same patient obtained on the same scanner. Intensity standardization was implemented in MATLAB R2014b (Mathworks, Natick, MA) using the method presented in Madabhushi et al. 28 , Another MRI artifact, bias-field inhomogeneity manifests as a smooth variation of signal intensity across the structural MRI, and has been shown to significantly affect computerized image analysis algorithms. Bias field artifacts were corrected for by means of the popular N4 bias-correction method 29 , which incrementally de-convolves smooth bias field estimates from acquired image data, resulting in a bias field corrected image. Similarly, hyperintense FLAIR signals correlate with greater interstitial leakage and low cellular density, reflecting edema. Therefore, T2w and FLAIR scans were used to identify edema and necrosis and enhancing tumour was delineated based on Gd-T1w MRI.

Segmentation.
To evaluate the effect of inter-observer variability in contouring the tumor sub-compartments, two experts were asked to independently manually segment 20 randomly chosen cases of GBM from our cohort (6 long term cases, 7 mid-term cases and 7 short term cases). Both the radiologists were provided separately with treatment naïve scans of Gd-T1w, T2w and FLAIR protocols. We obtained an average Dice Similarity Coefficient (DSC) scores across the 20 cases from the enhancing and edematous/nonenhancing region to be over 0.80. Compartment-specific radiomic feature extraction from MRI scans. A total of 30 2D radiomic features were extracted individually from every sub-compartment (edema/nonenhancing tumor, necrosis, enhancing tumor) for each of the 3 MR protocols (Gd-T1w, T2w, FLAIR). This resulted in a total of 270 features extracted for every study. The feature set for every study included 5 Laws energy, 12 Gabor, and 13 Haralick features on a per-pixel basis. A median feature value was then calculated from the feature responses of all pixels within the region of interest. All feature calculations were performed using in-house software implemented in MATLAB R2014b platform. A brief description of the extracted radiomic features is as follows: Detailed description of the set of features employed in this work and its possible relationship to the pathophysiology of GBM is provided in Table 2 below. The complete list of features extracted has been provided in the Supplementary spreadsheet (Sheet 1). processes (the hypoxia pathway in our case) and calculates an enrichment score for every patient in the cohort when paired with the 21 hypoxia associated genes 37 . Figure 2 shows the unsupervised clustering of the 21 gene set using euclidean distance, which identified three clusters low, medium, and high, using the HES values. The range of HES for each of the clusters is as follows: HES low ≤ 3141.99, 3142 ≤ HES Mid ≤ 4102.99, and HES High ≥ 4103. HES for each of the 85 GBM patients has been provided in the Supplementary spreadsheet (Sheet 2). Experimental Setup. Experiment 1: Identifying MRI surrogate markers of tumor hypoxia pathway. A total of 115 patients were split into 85 cases for training and 30 were held-out for validation. In the training phase, 85 patient studies were setup in a 3-fold stratified cross-validation for 50 iterations to create a total of 150 sets such that the samples in the training set were class-balanced. Each set consisted of randomized two-thirds data sampled into training and one-third data used for testing. For every training set, we obtained a mutual information score (captures linear or non-linear mutual dependence between two independent variables) between the radiomic features with the hypoxia enrichment score, and ranked the features based on their mutual information scores during each cross-validation run 21 . Subsequently, we obtained the frequency of occurrence of every feature across 150 runs of cross-validation (50 runs of 3-fold cross-validation). A subset of eight radiomic features with highest frequency of occurrences across the cross-validation runs were retained, for further analysis.    Statistical Analysis. Survival curves were compared statistically by a Cox proportional hazards model. All statistical analyses were performed using the survival package in R 39,40 . Hazard ratios (HR) were used to quantify the effect of individual feature on survival. Features yielding negative regression coefficients (i.e. low feature values correlated with long term survival) in our Cox Model produce a HR between 0 and 1; features yielding positive regression coefficients (i.e. low feature values correlated with short term survival) produce a HR between 1 and infinity. We also computed Concordance indices (C indices or C statistic) for each of our univariate and multivariate analysis experiments in R. C indices is the fraction of all pairs of subjects whose predicted survival times are correctly ordered (i.e. concordant with actual survival times). C indices = 1 indicates that the model has perfect predictive accuracy, and C indices = 0.5 indicates that the model is not better than random chance. Table 3 lists the top 8 radiomic features that were identified to be most associated with the hypoxia enrichment score using the mutual information feature selection method. The most associated radiomic features included Laws energy (R5R5, E5E5, S5S5) from enhancing tumor and edema capturing ripples, edges and spots. The top 8 features also included entropy, difference variance, and energy features from the Haralick family, which capture structural heterogeneity within the image texture, extracted from edema/nonenhancing and enhancing tumor. Figure 3 shows a single 2D slice of the original Gd-T1w MRI scan with annotations of the 3 tumor sub-compartments (necrosis outlined in green, enhancing tumor in yellow, and edema in brown), and the corresponding Haralick feature map for three different patients with high, medium and low hypoxia enrichment scores, respectively. While the original Gd-T1w MRI scans may not be able to visually capture the underlying tumor heterogeneity of hypoxic tumors, the radiomic features were found to be distinctly different across the varying degrees of hypoxia (low, medium, and high) (Fig. 3). High radiomic feature expressions corresponded to high hypoxia enrichment score (Fig. 3 (f)), while low feature expressions corresponded to low hypoxia enrichment score (Fig. 3 (d)).

Employing radiomic surrogate markers of tumor hypoxia in predicting patient survival in GBM.
The top 8 radiomic surrogate markers of hypoxia (as identified in Experiment 1) were also found to be significantly associated with survival (Fig. 4). Within the training set, Kaplan-Meier (KM) survival analysis between (a) STS versus LTS, and (b) MTS versus LTS showed a statistically significant separation between the survival curves as quantified via the log-rank test. Figure 4(a) shows the KM curve for patients with STS and LTS (p = 0.0056) while 4(b) shows the KM curve generated using the radiomic features between STS versus MTS (p = 0.8593), Fig. 4(c) shows the KM curves for MTS versus LTS (p = 9.2112×10 −6 ). On the validation set, significant differences in KM curves were observed for the short-, versus long-term survival patients (p = 0.0032) (Fig. 4 (d)  a C-index of 0.74. While, the KM curves for mid-term and long-term survival were not found to be statistically significantly different (p = 0.2093), the C-index was found to be 0.73. Table 4 lists the hazard ratio and concordance indices for the clinical parameters and combined radiomic features for distinguishing STS versus LTS, and MTS versus LTS. Interestingly, we found that combining clinical and radiomic features improved the concordance index in predicting overall survival in STS versus LTS (C-index = 0.69 and 0.83 on training and validation set respectively) and MTS versus LTS (C-index = 0.7 and 0.81 on training and validation set respectively) as compared to using either clinical features or radiomic features alone (Table 4).

Discussion
Currently, there is a lack of well validated non-invasive biomarkers that can predict the extent of hypoxia and potentially stratify patients that may be more suited for adjuvant anti-angiogenic treatments. For example, pre-clinical studies have shown that Bevacizumab (a humanized anti-VEGF monoclonal antibody for GBM therapy) attempts to normalize tumor vasculature, reduce hypoxia, and improve drug delivery 41 . Another challenge post-anti VEGF therapy is that patients sometimes fail treatment due to multiple intrinsic properties of the tumor (such as up-regulation of multiple pro-angiogenic factors, enhanced invasion and migration), where a GBM with high levels of hypoxia would indicate that the tumor has found other pathways to proliferate (for example, using ANG2) 7 . Therefore identifying non-invasive methods to monitor the hypoxic micro-environment has significant clinical implications in designing personalized treatment, as well as monitoring response to treatment in GBM patients. In this study, we investigated the relationship between MRI based radiomic features obtained from tumor sub-compartments and the corresponding gene expression data of the hypoxia pathway, as observed on the HES. We found prognostic radiomic features that were capable of distinguishing low, medium, and high levels of hypoxic extent (as defined from the HES), which could potentially also serve as imaging surrogates of overall survival in GBM.

Radiomic surrogate markers of hypoxia enrichment score. Our study identified Law energy and
Haralick features from the edematous/nonenhancing tumor and enhancing tumor region on FLAIR and Gd-T1w MR sequences to be highly associated with the hypoxia enrichment score. We believe that in the enhancing region of GBM, tortuous vessels of hypoxia induced neo-angiogenesis pile up and bulge to manifest as ripples and organization of pseudopalisades in the immediate vicinity of necrosis, evince as hypo-intense rings or spots and therefore are captured by law energy features on FLAIR 42 . This finding is concordant with Diehn et al., who found that their genomic expression module of hypoxia correlated with the enhancing region (p = 0.012) on Gd-T1w images 20 . Barajas et al., also found that enhancing regions of GBM with elevated relative cerebral blood volume (rCBV) were significantly correlated with hypoxia (p < 0.02) 43 . Similarly, we identified Haralick features (Entropy, difference variance) from the edematous region to be correlated with hypoxia. These Haralick features might potentially be quantifying the structural heterogeneity in highly hypoxic GBM tumors as observed on Gd-T1w MR sequences. For example, high values of entropy feature reflects high diversity in grey levels of diverse group of pixels (thus quantifying image heterogeneity) 12 .
Role of hypoxia-radiomic surrogate markers in predicting short-, medium-, long-term survivors of GBM. Imaging 12 . Similarly, in concordance with previous findings 12,47 , the combination of radiomic features with clinical parameters (Table 4) were found to improve prediction of GBM survival, as compared to radiomic features and clinical parameters alone. This suggests the prognostic potential that hypoxia-associated radiomic features offer in conjunction with clinical parameters in characterizing overall survival in GBMs. Lastly, it was observed that the short and medium term survivors were not separable on both the training and validation cohort (Fig. 4), suggesting that these categories may potentially represent a more aggressive radiomic phenotype compared to the long term survival cases.
The work presented in this study did have its limitations. In this study, we limited our analysis to only identifying associations of radiomic features with the hypoxia enrichment score, due to hypoxia's involvement in chemo-radiation resistance and poor outcomes in GBM. However, the radiomic features that were found to be associated with HES, may also be representative of other carcinogenic signaling pathways, such as tumor infiltration, proliferation, and angiogenesis that are implicated during hypoxia, and are known to contribute to poor outcome. For example, high expression of VEGFA (1 of the 21 genes that contributed to HES), promote repeated cycles of neo-angiogenesis that lead to microvascular hyperplasia, proliferation, and invasion in GBM tumors 48 . An extensive analysis of the association of the radiomic features with the other known pathways contributing to poor outcome (e.g. tumor infiltration, proliferation, and angiogenesis) will be a part of future study. Additionally, the results presented in this paper are preliminary and constrained by a relatively small sample size. A larger independent validation of the radiomic surrogate markers of hypoxia will need to be performed to further validate our preliminary findings. As the data was retrospectively collected from TCIA, another limitation of our study is that only routinely acquired MRI sequences (Gd-T1w, T2w and FLAIR) were used for analysis, and did not employ any advanced imaging (i.e. perfusion, DWI) including PET imaging.

Conclusion
In this study, we investigated the feasibility of computer-extracted radiomic features from different sub compartments of the tumor on treatment-naïve routine MRI in predicting extent of hypoxia and overall survival in GBM patients. The results suggest that radiomic features from the enhancing and edematous regions appear to be predictive of the extent of hypoxia (as observed using hypoxia enrichment score on mRNA data). The radiomic features on the validation set were also found to be prognostic of LTS vs STS. The identified radiomic features in this work could be used to monitor hypoxia, help determine timeline of treatment resistance, and evaluate the efficacy of anti-angiogenic therapy in GBM and other tumors.
While in this work, we limited the radio-genomic analysis to capturing radiomic imaging phenotypes that were associated with hypoxia pathway, our presented radiogenomic framework, could potentially in the future help identify imaging biomarkers associated with other key pathways such as tumor infiltration, proliferation, and angiogenesis, that are implicated in GBM, and contribute to poor outcomes. Future work will also focus on incorporating complementary imaging parameters obtained from advanced imaging (PET, perfusion, DWI) which may further improve survival prediction using radiomic analysis, while taking into account the extent of resection and subsequent treatment.