Multiparameter MRI Predictors of Long-Term Survival in Glioblastoma Multiforme.

Standard-of-care multiparameter magnetic resonance imaging (MRI) scans of the brain were used to objectively subdivide glioblastoma multiforme (GBM) tumors into regions that correspond to variations in blood flow, interstitial edema, and cellular density. We hypothesized that the distribution of these distinct tumor ecological "habitats" at the time of presentation will impact the course of the disease. We retrospectively analyzed initial MRI scans in 2 groups of patients diagnosed with GBM, a long-term survival group comprising subjects who survived >36 month postdiagnosis, and a short-term survival group comprising subjects who survived ≤19 month postdiagnosis. The single-institution discovery cohort contained 22 subjects in each group, while the multi-institution validation cohort contained 15 subjects per group. MRI voxel intensities were calibrated, and tumor voxels clustered on contrast-enhanced T1-weighted and fluid-attenuated inversion-recovery (FLAIR) images into 6 distinct "habitats" based on low- to medium- to high-contrast enhancement and low-high signal on FLAIR scans. Habitat 6 (high signal on calibrated contrast-enhanced T1-weighted and FLAIR sequences) comprised a significantly higher volume fraction of tumors in the long-term survival group (discovery cohort, 35% ± 6.5%; validation cohort, 34% ± 4.8%) compared with tumors in the short-term survival group (discovery cohort, 17% ± 4.5%, P < .03; validation cohort, 16 ± 4.0%, P < .007). Of the 6 distinct MRI-defined habitats, the fractional tumor volume of habitat 6 at diagnosis was significantly predictive of long- or short-term survival. We discuss a possible mechanistic basis for this association and implications for habitat-driven adaptive therapy of GBM.

Introduction

Glioblastoma multiforme (GBM) typically exhibits substantial intratumoral heterogeneity at both microscopic and radiological spatial scales (1). Analysis of genomic patterns from The Cancer Genome Atlas (TCGA) database led to a general molecular model that identified 4 distinct “species” of GBM: proneural, neural, classical, and mesenchymal (2). However, more recent studies (3) found substantial spatial variations, so that, in some cases, all 4 species could be observed in different regions of the same tumor. Canoll et al. used RNA-sequencing and histological analysis of image-guided biopsies to show differences in cellular and molecular markers between tissue taken from the contrast-enhancing (CE) core and that from the nonenhancing (NE) margins of GBM tumors (4). Characteristic metabolic differences between the CE and NE regions in GBM have also been identified by 1H magnetic resonance spectroscopy (5). Machine learning on patterns in standard brain magnetic resonance imaging (MRI) images, and parameter maps from diffusion tensor imaging, and dynamic susceptibility contrast-enhanced (DSC)-MRI have been reported to correlate with molecular subtype and survival in newly diagnosed patients with GBM (6). Radiogenomic analysis informed by spatially localized biopsies has identified spatially complex distributions of molecularly distinct subpopulations in GBMs (7). Although such spatial variations in expression of molecular and pathologic markers, metabolism, and radiologic imaging patterns are known to exist in all solid tumors, the origin and the clinical significance of this heterogeneity remain subjects of investigation.

Heterogeneity within tumors may drive resistance to both untargeted and targeted therapies (8). Reliance on conventional maximum tolerated dose–based treatment regimens may accelerate the unopposed proliferation of resistant populations by eliminating the susceptible populations and the attendant competition for space and substrate. Enriquez-Navas et al. recently showed that an evolution-based adaptive therapeutic strategy that exploits such competition between subpopulations of tumor cells could prolong progression-free survival in preclinical models of breast cancer (9). An ongoing clinical trial in prostate cancer (10) has shown that evolutionary dynamics can be successfully integrated into clinical cancer treatment protocols, and it highlighted the unmet need for noninvasive metrics of intratumoral subpopulation changes during treatment.

In the present work, we build upon a conceptual model of GBMs as spatially heterogeneous complex adaptive systems in which tumor growth and response to therapy are governed by eco-evolutionary interactions between the tumor microenvironment and phenotypic properties of local cellular populations. This model posits an explicit and predictable link between macroscopic tumor features observed radiologically and the molecular-, cellular-, and tissue-scale properties of the underlying cancer cell populations. In this model, we hypothesize that radiologically apparent spatial heterogeneity within each GBM can be quantified by some combination of a small number of distinct eco-evolutionary “habitats,” each of which may have different patterns of growth and invasion and may respond differently to therapy (11). Our approach builds upon methods developed in landscape ecology to bridge spatial scales. For example, field biologists are often tasked with estimating species distribution within a large area such as a county or state. Methods developed in landscape ecology typically begin with an analysis of satellite imagery of the region. By combining image channels containing nonoverlapping information (RADAR, infrared and visible light, for example), the biologist can divide the whole region into a patchwork collection of distinct habitats. By sampling the species distribution within each distinct habitat, the geographic distribution of each species over the entire region can be estimated (12, 13).

Multispectral clustering on MRI images has been used before to quantify spatial variations within tumors. Vannier et al. recognized the analogy between multispectral remote-sensing satellite imagery and multiparametric MRI and showed that signatures for “scene components” in the radiologic images could be computed (14–16). This approach can be used to further objectively subdivide the tumor itself into spatially distinct subregions (“habitats”) that harbor distinct subpopulations of tumor cells (11, 17, 18). Spatial heterogeneity of GBMs at radiological scales presents as regional variations in contrast enhancement and edema, and we have used multispectral clustering to decompose each glioma into a small number of distinct “habitats” based on their intensity on different MRI sequences. Tumor voxels were clustered by the calibrated signal intensities on contrast-enhanced T1-weighted (T1W-CE) and fluid-attenuated inversion-recovery (FLAIR) sequences into 6 distinct “habitats” based on low- to medium- to high-contrast enhancement and low–high signal on FLAIR scans. The long-term survival (LTS) cohort (>36 months postdiagnosis) were found to have a significantly higher fraction of habitat 6 (high CE and high FLAIR signal intensity) compared with the short-term survival (STS) cohort (≤ 19 months postdiagnosis) in both the discovery and validation cohorts. We discuss a possible mechanistic basis for this association between habitat 6 and survival in GBM, and implications for habitats-driven adaptive therapy of GBM.

Materials and Methods

Discovery Cohort

In this work, we have used the terms “discovery” (or training) and “validation” as they are understood in the field of machine learning, namely, to refer to the specific steps of training–validation–test in model development (19). Following IRB approval, patients with pathologically confirmed primary GBM and available preoperative T2-weighted (T2W), FLAIR, unenhanced T1W, and T1W-CE scans were identified retrospectively from a single participating institution. Median survival in glioblastoma is reported to be between 12 and 18 months postdiagnosis (20, 21). Recent estimates of 5-year survival rates for patients receiving maximal safe resection, concurrent radiotherapy and chemotherapy, and adjuvant chemotherapy are ∼10% (22). Our original intent was to investigate MRI habitats in high-grade gliomas from subjects who survived >5 years postdiagnosis. However, after application of the additional requirement that certain MRI scans be available at diagnosis, we had to downgrade this criterion to >3-year survival postdiagnosis of GBM so as to form cohorts with reasonable numbers of subjects. Thus, an LTS group comprising 22 subjects who survived >36 months postdiagnosis (median survival, 62.6 months; range, 36–107 months) was created. A control STS group of 22 subjects who survived <19 months postdiagnosis (median survival, 11.6 months; range, 2.5–19 month) was created to individually match to LTS subjects on age and calendar year of diagnosis.

Validation Cohort

Following IRB approval, patients with pathologically confirmed primary GBM and available preoperative T2W, FLAIR, T1W, and T1W-CE MRI scans were identified retrospectively from a multi-institutional database, matching on age and sex. The LTS group included 15 subjects who survived >36 months postdiagnosis (median survival, 86.6 months; range, 39–177 months), while the STS group included 15 subjects who survived ≤19 months postdiagnosis (median survival, 12.6 moths; range, 1.8–19 months).

Patient Population Statistics

Additional demographic and clinical covariates of relevance to this study are shown in Table 1.

Table 1.

Demographic and Clinical Characteristics of Patients in the Discovery and Validation Cohorts According to LTS and STS Status

Characteristics	LTS	STS
	Discovery Cohort
	(N = 22)	(N = 22)
Median Age (years)	50.5 (range: 22–74)	50.5 (range: 28–72)
Percent Male	59.1	63.6
Percent College Graduate[a]	45.5	23.8
Median KPS Score[a]	90%	80%
Median Year Diagnosed	2010	2011
Percent Completed Stupp Protocol[b]	37	0
Median Survival (Months)	67.7 (range: 36–126)	11.5 (range: 2.5–19)
	Validation Cohort
	(N = 15)	(N = 15)
Median Age (years)	50 (range: 23–68)	62 (range: 23–78)
Percent Male	67	60
Median Education (years)	Unknown	Unknown
Median KPS Score	90[c]	90[d]
Median Year Diagnosed	2009	2009
Percent Completed Stupp Protocol	66.7	26.7
Median Survival (months)	86.6 (range: 39–177)	12.6 (range: 1.8–19)

a 1 STS missing education; 3 LTS and 6 STS missing KPS score.

b As defined in PubMed PMID: 15758009. Results based on 20 LTS and 16 STS patients with complete information on receipt of the chemoradiation protocol. A total of 7 patients underwent biopsy as the only form of surgery (1 LTS and 6 STS).

c 10 missing values.

d 11 missing values.

Image Registration

For each patient, the FLAIR, T1W, and T1W-CE images were coregistered with the T2W images using in-house MATLAB (MathWorks, Natick, MA) software (top panel in Figure 1). As part of this process, the FLAIR, T1W, and T1W-CE images were resampled to match pixel dimensions and slice thicknesses with the reference T2W images. Spatial alignment was performed using a combination of rigid and affine geometrical transformations.

Figure 1.

Fluid-attenuated inversion-recovery (FLAIR), T1-weighted (T1W), and contrast-enhanced T1-weighted (T1W-CE) images were coregistered with and resampled to match voxel dimensions in the reference T2W scans (top panel). A contour was manually drawn to circumscribe the CE tumor in all applicable slices on postregistration T1W-CE images (middle panel). Normal white matter and cerebral spinal fluid (CSF) were automatically segmented (middle panel, details in Figure 2). Voxel intensities were calibrated against white matter (WM) and CSF to permit cluster analysis of voxels pooled across patients on each type of magnetic resonance imaging (MRI) scan (middle panel). Pooled voxels from within the CE tumor mask were clustered into 6 habitats using the criteria listed in Table 1 (bottom panel). Also shown in the bottom panel is a 3D stack of maps of habitats 1–6 in an example tumor, for illustrative purposes.

Figure 2.

Automatic segmentation procedure to locate WM and CSF volumes within the brain for use in intensity calibration. Intensities within the T1W-CE − T1W difference volume (ΔT1W) of a given subject were clustered into low- and high-intensity classes by Otsu thresholding. A mask of voxels in the low-intensity class was applied to the T2W image and further subdivided into low- and high-intensity clusters by Otsu thresholding. The resulting mask of voxels in the low-intensity cluster was applied to the T1W image, which was again subdivided into low- and high-intensity clusters with the high-intensity class labeled as “normal white matter” (reference region 1). The mask of high-intensity voxels from the T2W image was applied to the FLAIR image, and it was again subdivided into low- and high-intensity clusters with the resulting low-intensity cluster labeled as “CSF” (reference region 2).

Tumor Segmentation

In this work, we restricted our analysis of intratumoral “habitats” to the CE portion of the tumor volumes. For this purpose, a contour was manually drawn to circumscribe the CE tumor in all applicable slices on postregistration T1W-CE images (middle panel in Figure 1).

Intensity Calibration

The next step in our image processing pipeline was intensity calibration (middle panel of Figure 1), the objective of which is to allow comparison of voxel intensities across patients on each given type of MRI scan. For this purpose, 2 reference normal tissue regions were automatically segmented as shown in Figure 2. In brief, intensities within the T1W-CE − T1W difference volume (ΔT1W) were clustered into low- and high-intensity classes using Otsu thresholding (23). Then, on T2W, voxels from the low-intensity class were subdivided further into low- and high-intensity clusters using Otsu thresholding. Voxels from the low cluster formed a volume of interest (VOI) that was applied to T1W, which was subdivided into low- and high-intensity clusters by Otsu thresholding, with the resulting voxels in the high-intensity class labeled as “normal white matter” (reference region 1). Voxels from the high T2W cluster formed a VOI mask that was applied to the FLAIR scan, and these were again subdivided into low- and high-intensity clusters using Otsu thresholding, and the low-intensity cluster was labeled as “CSF” (reference region 2). Voxel intensities on T2W, FLAIR, and precontrast T1W images were then linearly calibrated using “normal white matter” and “CSF” as reference tissues. The reference intensity values for these 2 tissues, respectively, were 81 and 183 on T2W, 587 and 464 on FLAIR, and 1099 and 748 on precontrast T1W, all in arbitrary units. These reference values were taken from the T2W, FLAIR, and precontrast T1W images of a patient chosen randomly from the discovery cohort, and do not carry any particular physiological meaning as such. Intensity calibration for T1W-CE was performed using the same linear transformation as computed for the associated precontrast T1W. Our input data comprise standard-of-care MRI images that were acquired with varying protocols per subject. Acquisition parameters such as the repetition time, echo time, and flip angle were not the same across all subjects for each scan type (T2W, FLAIR, T1W). Because MRI signal intensity is a nonlinear function of these acquisition parameters, linear calibration against 2 reference tissues may not necessarily be adequate for standardization of intensities per scan type. Fortunately, the range of excursions in these acquisition parameters across subjects was relatively small, and signal equation simulations indicated that calibration of raw signal intensity against 2 dissimilar reference tissues would provide satisfactory intensity calibration for other tissues with T1 and T2 values similar to or in-between those of the 2 reference tissues. The coefficient of variation of normal gray matter intensity across all patients was significantly smaller postcalibration as compared with precalibration on each of FLAIR, T1W, and T1W-CE images, and we took this to be evidence of successful intensity calibration (see online Supplemental Figure 1).

Multispectral Clustering to Define Intratumoral Habitats

Calibration of intensities per MRI scan type allows us to pool voxels over multiple patients for combined cluster analysis. This series of steps is depicted in the bottom panel of Figure 1. In brief, the manually drawn CE tumor mask was applied to the calibrated ΔT1W difference volume of each patient in the discovery cohort, and the voxels within the mask were pooled over all subjects and clustered by Otsu thresholding into 3 levels of contrast enhancement: CE1 (low enhancement), CE2 (medium enhancement), and CE3 (high enhancement). The low-, medium- and high-contrast enhancement thresholds identified on the discovery cohort were refined on validation, specifically that the maximum value of ΔT1W difference intensity was capped at 5000 arbitrary units postcalibration before Otsu thresholding. This was done to manage the skewing of the clustering process by a long 1-sided tail on the ΔT1W difference intensity histogram in some patients. Each of these 3 clusters was further subclustered into 2 classes around a calibrated value of 600 on FLAIR, a threshold value that is similar to the mean intensity of normal white matter over all subjects after calibration. The final habitat definitions are listed in Table 2.

Table 2.

Intratumoral Habitats' Definitions on Calibrated FLAIR and ΔT1W Intensities

	Calibrated FLAIR Image Intensity	Calibrated ΔT1W Difference Intensity
Habitat 1	≤600	≤303
Habitat 2	>600	≤303
Habitat 3	≤600	303 < ΔT1W ≤ 790
Habitat 4	>600	303 < ΔT1W ≤ 790
Habitat 5	≤600	>790
Habitat 6	>600	>790

Statistics and Survival Analyses

Absolute tumor volume, habitat volumes, and habitat volume fractions for each habitat were computed. Statistical analyses were performed using GraphPad Prism 7 (GraphPad Software, La Jolla, CA). Data normality was assessed using the D'Agostino–Pearson test, and significance of differences in habitat volumes between groups was assessed by 2-tailed unpaired t-tests. Survival analyses were performed using Kaplan–Meier survival curves, and statistical significance was computed using the log-rank test. For the Kaplan–Meier analysis, habitat volumes were dichotomized into 2 groups using the median score value.

Results

Mean tumor volumes at diagnosis were comparable between the LTS and STS groups in the discovery cohort (33 ± 6.6 cm3 vs. 37 ± 6.1 cm3, P = .62) (see online Supplemental Figure 2A). There was no statistically significant difference in mean tumor volumes at diagnosis between the LTS and STS groups in the validation cohort (33 ± 7.0 cm vs. 17 ± 4.8 cm3, P = .075), although there was a trend toward smaller tumor volumes in the STS group (see online Supplemental Figure 2B).

Figure 3 depicts differences in habitat 6 (high contrast enhancement and high FLAIR) content between a representative LTS subject (left; overall survival, 41+ months) and STS subject (right; overall survival, 3 months) at the time of tumor presentation before surgical intervention. In the discovery cohort habitat 6 comprised a significantly higher volume fraction (P < .03) of the tumor volume at diagnosis in long-term survivors (mean ± S.E.M. = 35% ± 6.5%; n = 22) compared with short-term survivors (mean ± S.E.M. = 17% ± 4.5%; n = 22) (Figure 4A). This finding was replicated in the validation cohort (P < .007), with habitat 6 comprising 34% ± 4.8% (n = 15) of the tumor volume in LTS subjects compared with 16% ± 4.0% (n = 15) of the tumor volume in STS subjects (Figure 4B).

Figure 3.

Habitat 6 (high enhancement and high FLAIR) on preoperative MRI comprises 23% of the tumor by volume in a long-term survivor (left, overall survival 41+ months) and 9% of the tumor by volume in a short-term survivor (right, overall survival 3 months).

Figure 4.

Habitat 6 (high enhancement and high FLAIR) was significantly higher in the LTS group relative to the STS group in both the (A) discovery and (B) validation cohorts.

Habitat 2 (low enhancement and high FLAIR) comprised a significantly lower volume fraction (P = .0126) of the tumor at diagnosis in long-term survivors (mean ± S.E.M. = 28 ± 5.7; n = 22) relative to short-term survivors (mean ± S.E.M. = 51 ± 6.8; n = 22) in the discovery cohort (Figure 5A), but this was not replicated in the validation cohort (Figure 5B). In parallel, habitat 1 (low enhancement and low FLAIR) was not found to be significantly different between LTS and STS subjects in the discovery cohort (Figure 5C) but comprised a significantly lower volume fraction (P = .0279) of the tumor at diagnosis in long-term survivors (mean ± S.E.M. = 3.2 ± 0.96; n = 15) relative to short-term survivors (mean ± S.E.M. = 12 ± 3.4; n = 15) in the validation cohort (Figure 5D). Minor inconsistencies in FLAIR intensity calibration across the patients may be the root cause of this variable finding, given that Habitats 1 and 2 belong to the low and high FLAIR clusters, respectively.

Figure 5.

Habitat 2 (low enhancement and high FLAIR) was significantly lower in the LTS group relative to the STS group in the discovery cohort (panel A), but this difference was not recapitulated in the validation cohort (panel B). Habitat 1 (low enhancement and low FLAIR) was not significantly different between the LTS and STS groups in the discovery cohort (panel C), but was significantly lower in the LTS group in the validation cohort (panel D).

Habitat 3 (medium enhancement and low FLAIR), habitat 4 (medium enhancement and high FLAIR), and habitat 5 (high enhancement and low FLAIR) were not significantly different between the LTS and STS groups in either the discovery or validation cohorts (see online Supplemental Figure 3).

Median percent of tumor volume occupied by habitat 6 in the discovery cohort (5.77%) was used as a cutpoint to dichotomize patients into high and low habitat 6 fraction groups. Kaplan–Meier survival analyses were then carried out separately in the discovery and validation cohorts using the prespecified cutpoint (5.77%) established for the discovery cohort, as a stringent test of reproducibility. Based on the median cutpoint, low and high fractions of habitat 6 were not associated with overall survival in the discovery cohort (Figure 6A; P = .62), but were statistically significant with respect to overall survival in the validation cohort (Figure 6B; P = .0001). In the discovery cohort, Kaplan–Meier 3-year survival rates were 45% and 55% in the < median versus ≥ median subgroups, respectively. In the validation cohort corresponding 3-year survival rates were 18% and 68%, respectively.

Figure 6.

Kaplan–Meier plots of overall survival in the discovery cohort (A). Survival of patients with habitat 6 volume fraction ≥ median (5.77%, n = 22) and < median (n = 22). Kaplan–Meier plots of overall survival in the validation cohort (B). Survival of patients with habitat 6 volume fraction ≥ median from the discovery cohort (5.77%, n = 19) and < median (n = 11).

Discussion

The overall goal of our work is to develop noninvasive imaging biomarkers that can be used to drive evolution-based adaptive therapeutic strategies for GBM. For any biomarker to be clinically useful, it must be computable reliably and reproducibly (24). MRI parameters such as ADC, T1, and T2, and with some limitations, also model-dependent parameters such as relative cerebral blood volume (rCBV), relative cerebral blood flow, and Ktrans, are comparable between data sets when standardized protocols are utilized (25–35). Parameter maps are therefore attractive for computing tumor habitats consistently across patients and scan dates, but these maps are not routinely collected as part of standard-of-care imaging. The subjects in our study received their initial diagnostic scans at a variety of institutions including at community radiology facilities, as a result of which there was great variability in the type and quality of scans that were available for retrospective analysis. In particular, we were unable to curate sufficient numbers of LTS subjects with available ADC maps at diagnosis. We therefore sought to compute intratumoral habitats using FLAIR, T1W, and T1W-CE scans after calibrating raw MRI pixel intensities against 2 reference tissues.

High signal on ΔT1W is indicative of either good perfusion or high microvascular leakiness. High intensity on FLAIR images in glioma represents a mixture of vasogenic edema, which arises from leakage of plasma into regions with low cell density, and tumor cell infiltration along long white matter tracts (36). Our retrospective study shows, in both a discovery cohort and a validation cohort, that tumors in LTS subjects have a significantly higher fraction of habitat 6 (high contrast enhancement and high FLAIR signal intensity) than STS. Particularly striking is the similarity in habitat 6 content of LTS tumors between the discovery and validation cohorts (35% and 34%, respectively) and of STS tumors between the discovery and validation cohorts (17% and 16%, respectively). We divided tumor regions with high signal intensity on ΔT1W calibrated difference images into 2 distinct habitats with either high or low FLAIR signal. Low FLAIR signal would be expected in regions with high contrast enhancement stemming from good perfusion, which would be conducive to high cellular density, although not necessarily where the enhancement arises from microvascular leakiness. Our results demonstrate the high contrast enhancement and high FLAIR signal habitat is strongly associated with patient survival.

In a preliminary study of pretreatment MRI examinations from 32 patients with GBM enrolled in the TCGA, Gatenby et al. showed that GBM tumor habitats defined on FLAIR and T1W-CE images could be used to differentiate patients who survived <400 days from patients who survived ≥400 days postdiagnosis (37). A follow-up study indicated that incorporating information from 3 MRI sequences, namely, T2W, FLAIR, and T1W-CE, improved prediction of survival time in patients with GBM (38). LaViolette et al. similarly clustered voxels into low, medium, and high classes on T1W, T1W-CE, FLAIR and apparent diffusion coefficient of water (ADC) maps to divide GBM tumors into 81 habitats, and identified 5 specific habitats that when present at higher volumes correlated with poorer prognosis (39). Recently, Juan-Albarracín et al. analyzed preoperative DSC-MRI and FLAIR scans of 50 patients with GBM to compute tumor habitats on the basis of rCBV, relative cerebral blood flow, and edema, and they found a surprising correlation between longer survival times and lower indices of perfusion (40). Boonzaier et al. report that tumor habitats reflecting low ADC values intersecting with high rCBV values demonstrate a significantly elevated choline-to-N-acetylaspartate ratio on 1H magnetic resonance spectroscopy, and that a higher proportion of this habitat within the NE region of GBM is associated with poor overall survival (41). Interpatient diversity in overall imaging patterns of growth and invasion has been associated with tumor aggressiveness and clinical outcomes across patients (42–45). Our investigation leverages unique resources of data including patients with exceptionally long follow-up for prognosis in glioblastoma.

Standard-of-care therapy in newly diagnosed GBM is maximal safe surgical resection followed by concomitant radiation therapy and temozolomide for 6 weeks, followed by adjuvant temozolomide for 6 monthly cycles (46). Thereafter, subjects in our retrospective study would each also have received a variety of investigational and/or palliative treatments, including extended cycles of temozolomide. Our findings suggest that one or more characteristics of the radiologically visible initial tumor mass define an intrinsic prognostically relevant tumor feature that continues to influence patient outcome months, and even years, after diagnosis. It is possible that the radiologic appearance of habitat 6 is a shared feature of disparate favorable markers in GBM, such as Isocitrate Dehydrogenase (IDH) mutation status (47), mesenchymal subtype (48) or lymphocyte cytokines such as CXCR4 (49). Alternately, one can hypothesize that components of the immune system in the LTS subjects retain the ability to recognize tumor antigens present in the original mass that are retained in the recurrent mass. Immune infiltrates in the tumor would be consistent with the MRI characteristics of habitat 6, namely, high contrast-enhancement and high tumor-associated edema. Pathological studies have shown that increased CD8+ T cell infiltrates in newly diagnosed GBM is associated with long-term survival (50), and we hypothesize that increased FLAIR signal in well-perfused—and presumably cellular—regions may be indicative of interstitial edema related to inflammatory changes caused by an immune response. A definitive biological interpretation of our finding requires further investigation.

Known weaknesses in our study include that the numbers in each survival group stratum were small and statistical power correspondingly limited to detect all but strong associations in the data. Specifically, while our analysis detected a significant difference between the LTS and STS groups in both the discovery and validation cohorts (Figure 4), on an individual patient basis, we could observe survival differences by a binary analysis around the median habitat 6 content in only the validation cohort (Figure 6). The need to improve calibration of raw MRI image intensities is revealed in the inconsistent significances of Habitats 1 and 2 in the discovery and validation cohorts (Figure 5). Additional covariates may also impinge upon our analysis. For example, in the discovery cohort, LTS and STS subjects were matched for parameters such as patient age and year of diagnosis, but LTS patients were nonetheless more educated and more likely to survive the completion of standard treatment. In the validation cohort, the LTS and STS groups were not matched for patient age and treatment regimens. It is unclear how these group differences might explain the present findings.

Only about 5% of patients with GBM undergoing standard of care survive ≥5 years postdiagnosis (46). Investigation of a cohort of rare long-term survivors identifies a “habitat” on initial multiparametric MRI scans that is significantly different than in a control cohort. Our working hypothesis is that habitat 6 corresponds to a microenvironment that selects for glioma cells that are either innately less aggressive or are more amenable to control by tumor-infiltrating leukocytes. Habitat imaging has the potential to provide noninvasive longitudinal biomarkers of intratumoral evolutionary and ecological dynamics for the informed application of adaptive therapy to manage GBM.