Higher vascularity at infiltrated peripheral edema differentiates proneural glioblastoma subtype

Background and purpose Genetic classifications are crucial for understanding the heterogeneity of glioblastoma. Recently, perfusion MRI techniques have demonstrated associations molecular alterations. In this work, we investigated whether perfusion markers within infiltrated peripheral edema were associated with proneural, mesenchymal, classical and neural subtypes. Materials and methods ONCOhabitats open web services were used to obtain the cerebral blood volume at the infiltrated peripheral edema for MRI studies of 50 glioblastoma patients from The Cancer Imaging Archive: TCGA-GBM. ANOVA and Kruskal-Wallis tests were carried out in order to assess the association between vascular features and the Verhaak subtypes. For assessing specific differences, Mann-Whitney U-test was conducted. Finally, the association of overall survival with molecular and vascular features was assessed using univariate and multivariate Cox models. Results ANOVA and Kruskal-Wallis tests for the maximum cerebral blood volume at the infiltrated peripheral edema between the four subclasses yielded false discovery rate corrected p-values of <0.001 and 0.02, respectively. This vascular feature was significantly higher (p = 0.0043) in proneural patients compared to the rest of the subtypes while conducting Mann-Whitney U-test. The multivariate Cox model pointed to redundant information provided by vascular features at the peripheral edema and proneural subtype when analyzing overall survival. Conclusions Higher relative cerebral blood volume at infiltrated peripheral edema is associated with proneural glioblastoma subtype suggesting underlying vascular behavior related to molecular composition in that area.


Introduction
In the late years, Central Nervous System tumor classification has shifted from being based on microscopic similarities between cells and their levels of differentiation [1] to additionally include genetic-based features [2]. This is particularly the case for glioblastoma, where several classifications have been defined: on the one hand, the World Health Organization (WHO) classification which distinguishes between IDH-wildtype and IDH-mutant glioblastomas [2][3][4] and, on the other, the Verhaak classification [5], consisting of 4 subtypes depending on mutations and molecular profile of various cancer-related genes. These subtypes are the mesenchymal, classical, neural and proneural, the latter being related to IDH mutations [5,6]. These new classification paradigms have improved the estimation of prognosis [7,8] and proposed specific therapeutic targets [9][10][11][12], especially for patients with proneural and mesenchymal type glioblastoma.
Considering that Magnetic Resonance Imaging (MRI) perfusion biomarkers have been associated with patients' overall survival [13][14][15] and cellular features [16,17], several studies were performed to analyze if there was a relationship between vascular biomarkers and the genomic subtypes classifications. Barajas et al. studied the influence of glioblastoma genetic and cellular features over MRI, concluding that they could spot the most malignant regions within the tumor [18]. Jain et al. demonstrated that combining Verhaak subtypes with vascularity markers at the enhancing tumor provides additional information as a survival predictor [19]. However, they found that the enhancing and non-enhancing regions of the tumor did not present any significant correlations with the genomic subclassification. Another study proposed that tumor blood volume determined by dynamic susceptibility contrast MR perfusion imaging was related to EFGR and to PTEN expression in some patients [20].
Most of these studies focus on the vascularity of the enhancing tumor region and only a few remarked the influence of the non-enhancing part of the tumor including the edema region [18,20]. Gill et al. [21] found molecular differences between Verhaak subtypes performing MRI-localized biopsies in the peripheral edema region. Similarly, Price et al. [22] discovered that metabolic and perfusion changes in this region could be found using multimodal MR images. In this sense, we hypothesize that the vascular parameters in the invasive margins of glioblastoma could be related to characteristic combinations of mutations.
The purpose of this article is to assess the correlation between the vascularity present at the infiltrated peripheral edema habitat at preoperative stage and Verhaak molecular classification. To do so, we propose the use of a multicentrically validated [23] automatic open service named ONCOhabitats (https://www.oncohabitats.upv.es) proposed by Juan-Albarracín et al. [15,24,25]. To ensure the comparability of our study, the analysis was performed on the TCGA-GBM open database [26], which contains MR images and molecular information. In the end, we found correlation between peripheral edema vascularity and specially the proneural glioblastoma subtype.

Patient selection
Our study included retrospective patients with glioblastoma from The Cancer Imaging Archive-TCGA-GBM [26], an open dataset which only uses anonymized data, available at: https://wiki.cancerimagingarchive.net/display/Public/TCGA-GBM (Version 3). The database consists of 262 histopathological validated glioblastoma patients, 66 of which had preoperative dynamic susceptibility contrast enhanced T2 � -weighted perfusion (DSC) imaging information. Three of them were excluded because they did not have genomic information available.
The remaining 63 belong to two different institutions with the following distribution: 48 in the first and the rest in the second. From the first institution, 6 were excluded because of poor perfusion acquisition, mainly due to having an incomplete field of view in the DSC images. Additionally, 5 were excluded due to post-processing errors when performing DSC quantification. From the second one, only 2 were not considered for having an incomplete FOV.
According to the Verhaak molecular classification [27,28], the group of patients would be divided into 10 classical, 17 mesenchymal, 11 neural and 12 proneural subtype glioblastomas, attending to the mutations and markers they presented. The cohort clinical and genetic data along with the subtype of each subject can be found in the S1 Table, whose complete information is retrieved from the original at the TCGA-GBM website [26].

MRI imaging acquisition
From both institutions only 17 studies were obtained using 3T magnetic resonance imaging machines, all belonging to the first institution. For the rest of them, 1.5T imagers were used.
DSC perfusion MRI was performed during the injection of the gadolinium-based contrast (0.1 mmol/kg) using 95 dynamics for the first institution and 60 dynamics for the second institution of T2 � -weighted gradient echo echoplanar images. T1, T2 and fluid attenuation inversion recovery (FLAIR) MR images were acquired before the DSC perfusion. Gadolinium enhanced T1-weighted (T1-Gd) images were obtained after. The acquisition matrices were either 128x128 or 256x256, with slice thicknesses ranging from 1.5 to 6.0 mm. S2 Table shows repetition time (ms) / echo time (ms) / flip angle (˚) of the acquisitions for every patient.

Computing vascular habitats
All the cases were processed using the Hemodynamic Tissue Signature service found in the ONCOhabitats platform [24]. It provides a reproducible [23] and automated methodology to define habitats withing the glioblastoma, based on the vascular properties of the lesion, as shown in Fig 1. Morphological sequences are used for segmentation of enhancing tumor and edema. The resulting segmentation together with DSC perfusion maps are used to obtain the vascular habitats.
MRI preprocessing includes several modules for denoising, inhomogeneity correction, registration, brain extraction, motion correction or intensity standardization.
DSC perfusion maps are calculated following standard methods published in the literature. The block-circulant SVD decomposition method and numerical integration of gamma-variate concentration-time curves are used to estimate rCBF and rCBV respectively. Boxerman leakage correction method [29] is also employed to correct for T1-and T2-extravassation effects.
Glioblastoma tissue segmentation is performed through a Convolutional Neural Network (CNN) trained with more than 300 cases manually segmented by several expert radiologists. The CNN works with 3D MRI patches of post-Gadolinium T1-weighted, T2-weighted and FLAIR images. The network follows a U-net architecture with skip-connections and residualinception blocks, achieving a Dice performance in whole tumor segmentation of 0.90 [30]. Vascular habitats are obtained through an unsupervised methodology based on NL-SVFMM [31] method. A mixture of 4 classes is fit into the region delineated by the morphological segmentation of the enhancing tumor and edema to obtain the habitats. Several physiological constraints are introduced during the clustering process to drive the segmentation to a physically plausible delineation of the habitats. For the case of the IPE habitat, which is the relevant habitat studied in this work, their apparition is constrained, following the specifications of Guo et al. [32], to a 2 cm margin around the enhancing tumor, which is considered the plausible region where edema could be infiltrated by the tumor. The remaining habitats, that is, the high angiogenic tumor (HAT), low angiogenic tumor (LAT) and the vasogenic peripheral edema (VPE), also introduce constraints limiting their locations and proportions of occurrence within the lesion. rCBV max , rCBV mean and rCBV median were calculated for each vascular habitat. rCBV max was defined as the 95 th percentile of the distribution of rCBV values within the ROI in order to increase robustness. Values of rCBV max , rCBV mean and rCBV median at each habitat for each subject are presented in the S3 Table.

Statistical analysis
Firstly, ANOVA and Kruskal-Wallis H test were performed at the ET -in order to have comparable results with the current literature-, consisting of HAT and LAT regions, and at the IPE using all perfusion parameters (rCBV max , rCBV mean and rCBV median ) across Verhaak subclasses. This will serve as a first approximation for establishing any significant divergences in vascularity values regarding the Verhaak subtypes. For deepening the analysis on the specific differences between the four subclasses, Wilcoxon-Mann-Whitney tests were executed in each habitat. The comparison was made between classical, mesenchymal, neural and proneural classes individually one against the other, and comparing each one of them against the three remaining. For significant experiments, ROC curves were drawn for threshold optimization.
These statistical tests were considered significant when p-values were under 0.05. To correct for multiple testing, Benjamini-Hochberg false discovery rate correction was carried out for every study. Analyses were carried out on a personal computer with MATLAB R2018a (Natick, Massachusetts, USA).
Finally, survival Cox proportional hazards analysis were carried out in order to assess the effect of Verhaak subtypes on Cox overall survival models based only on perfusion parameters. To this end, univariate survival models were fitted using only rCBV max at IPE and at ET. Then, univariate models using only each subtype were fitted. Finally, multivariate models with both rCBV max and each subtype as cofactor were studied. The consequences of subtype addition can be either worsening or improving the fitting, indicating that the subtype provides redundant information or not, that is, there is an association between subtype and vascularity at IPE or not.
For each Cox model, Hazard Ratio (HR) with 95% confidence intervals (CI95), area under ROC curve (AUC) and p-values are reported. Significance will be considered when p-values are under 0.05. Analyses were carried out on a personal computer using R statistical analysis software [33].

Verhaak subtypes and rCBV in vascular habitats
In Fig 2 the Box-Whiskers representation of rCBV max values at the vascular habitats for each Verhaak subtype is represented. Firstly, as expected, a decrease in vascularity can be observed as we move from HAT to VPE. The proneural subtype shows higher vascularity in every habitat. However, there is important overlap at the enhancing areas, whereas at IPE and, to a lesser degree, at VPE the difference with the rest of the subtypes is bigger. This may point to vascularity differences in the peripheral region.
These results are consistent with the ones presented in Table 1. rCBV values at the ET were not significantly different among Verhaak subclasses neither performing an ANOVA nor a Kruskal-Wallis test. However, performing the same analyses on the IPE, significance is found for every rCBV metric (i.e. Max, Mean, Median). Results of subtype-specific tests at IPE are presented in Table 2. Comparing rCBV max values at the IPE habitat between the Verhaak subtypes, we obtained that the proneural tumor subtype has a significantly differentiated peritumoral vascularity. Note that the most significant difference was found between proneural and mesenchymal subtype. Mesenchymal subtype vascularity was significantly different from the other three subtypes only before multiple-test correction. Classical and neural subtypes had indistinguishable vascularity at the IPE according to this test, both one from the other and the two from the other subtypes. The same tests were performed for the rest of habitats in the S1 Appendix, yielding only significant results for proneural against mesenchymal rCBV max at VPE.

Proneural subtype differences in rCBV at IPE region
When comparing proneural rCBV mean and rCBV median at IPE against the rest of the subtypes, significant differences were also found (corrected p-values of 0.0428 and 0.0420 respectively). S1 Fig shows ROC curves for threshold optimization of the rCBV max at IPE using the three significant experiments after multiple test correction: (1) differentiating proneural from mesenchymal, (2) proneural from neural and (3) proneural from the of subtypes together. Optimal  value for distinguishing proneural from mesenchymal was a rCBV max of 3.10 at IPE, whereas for proneural from neural the best cutoff was a rCBV max of 3.01 at IPE. Finally, the optimal threshold for differentiating proneural from the rest of Verhaak subtypes was a rCBV max of 3.12 at IPE. Table 3 shows Cox proportional hazards regression for rCBV max at IPE and Verhaak subtypes. Vascularity at IPE alone is significatively associated with overall survival. Proneural subtype  also yields significant results. When adding subtypes as cofactors, HR, p-values and AUC are maintained relatively stable for every subtype except proneural. In the latter, p-values escalated to non-significance and HR decreased, affecting also confidence intervals. AUC increase is no longer meaningful as neither regressor is significantly associated with survival. All of this points to a blurring in the effect of the biomarker due to the addition of correlated redundant molecular information. The same test was performed for rCBV max at ET in the S2 Appendix. As a well-known biomarker, it also shows significant association with overall survival by its own. In this case, however, when adding the subtypes as cofactors, HR and significance is maintained for every subtype, including proneural, sometimes even improving AUC.

Discussion
In this study, we examined whether the vascular properties of the potentially infiltrated peripheral edema habitat were correlated with Verhaak molecular subclasses, and especially the proneural subtype. The results show that a value of rCBV max higher than 3.12 at IPE is significantly related to proneural subtype.
There are several studies aiming to determine the influence of genetic expression patterns in MRI features within the edema region. Carrillo et al. pointed out that edema could have prognostic importance in cases when MGMT promoter is methylated [34], which is a common trait in the proneural subtype [35]. Naeini et al. discovered that T2 and FLAIR volume hyperintensity representing edema was higher in proneural phenotypes [36] and Zinn et al. published that, by stratifying into high and low FLAIR radiophenotypes, they could identify glioblastoma subtypes [37]. Finally, the study of MRI perfusion and genetics of GBM from Barajas et al., pointed out the need for a deeper understanding of peritumoral non-enhancing tumor for its risk in future progression as its genetic expression pattern differs from that of the enhancing lesion [18].
In light of these results, in this study we analyzed the radiomic relevance of the edema using a more detailed characterization of edema heterogeneity by differentiating between nonenhancing (i.e. IPE) and vasogenic edema (i.e. VPE), based on ONCOhabitats approach [15]. This allowed to overcome a limitation pointed out in previous studies [36] and thus identify the IPE as a region with a radiomic relevance when studying the proneural type. In particular, we found a significant association of the rCBV at the IPE with the glioblastoma molecular profile. The prognosis potential of the rCBV in glioblastoma at the ET ROI has been extensively studied [14,38,39]. The prognosis potential of the rCBV in glioblastoma in the non-enhancing region was also suggested in some studies [15,23,40]. Jain et al. performed a survival analysis estimating HR for rCBV in both the enhancing and non-enhancing areas adjusting them for the Verhaak molecular subtypes and using the same database as this study [19]. They found that statistical significance for vascularity enhancing areas improved, suggesting additional information provided by molecular profile. However, for the edema region, it did not improve when adding the molecular information. The different molecular behavior of enhancing and non-enhancing regions is consistent with our results. Moreover, our method allowed to find bigger differences when adjusting for the proneural subtype. Thus, the Cox regression model performed in this study suggested that the predictive power of rCBV at non-enhancing areas could be related to its relationship with survival-related mutations.
As Verhaak noted when describing the molecular subclasses, they can be therapeutically relevant [5]. The fact that the biggest differences were found between the proneural and mesenchymal subtypes may point out that their vascular behavior in peritumoral regions varies broadly. As these two have been the most clinically relevant [10,12], being able to identify them by specific perfusion features can be crucial for diagnosing and treating glioblastoma patients.
Additionally, in 2016 the WHO published a glioblastoma classification mainly differentiating whether the IDH gene presents itself as mutant or wild-type [2]. Our findings could also agree with the classification: the proneural subtype was the most significantly different subtype from the three remaining and it has proven to be the most closely related to IDH1 mutations [5]. Unfortunately, not enough information was available to confirm if IDH-mutant vascularity at the peripheral areas was significantly different from wild-type.
Finally, in 2012 Strum et al. showed that the perceived better outcome for proneural subtype was caused by the presence of IDH-mutant patients and, when they were removed, proneural patients actually had the worst prognosis [41]. This is also supported by the studies that evaluated the TCGA-GBM dataset and found proneural had the shortest median survival [5,19]. Our findings show that the proneural subtype has the highest vascularity in the peripheral edema, which could be considered a surrogate for aggressiveness.
There were some limitations in this study. Firstly, though we were able to correlate the IPE habitat and the Verhaak subclasses, due to the retrospective nature of the study, we were not able to standardize MRI acquisition protocols. Secondly, it can be difficult to correctly assess the exact location of the infiltrated edema by a noninvasive manner. These could affect the potential relationships with molecular markers. Nonetheless determining the IPE habitat by an automated method for calculating the maximum of the cerebral blood volume to perform seems robust enough for our purpose, as shown by Á lvarez-Torres et al. [23]. Finally, despite having significant results, the sample size available in the dataset may be a statistical limitation.
Our study relies on the use of an automatic procedure to determine a more precise peritumoral ROI based on an open service [24] and the results can be replicated using the TCGA-GBM open dataset [26].

Conclusions
In conclusion, high IPE vascularity features are associated with the proneural subtype. Global vascularity differences between the four subtypes exist in this region especially due to proneural and mesenchymal influence. rCBV max at IPE is related to overall survival and carries specific molecular information.
Supporting information S1 Table. Clinical data of the final cohort. All clinical data is obtained from the TCGA-GBM open database [26]. Subtype and genetic information were retrieved from the UCSC Xena platform compilation [28] which is based on Brennan et al [27] (Supplementary Material 9). (XLSX) S2 Table. MRI scanners and protocols used for each subject. (XLSX) S3 Table. rCBV max , rCBV mean and rCBV median at every habitat for each subject.