Purpose:

Heterogeneity in glioblastomas is associated with poorer outcomes, and physiologic heterogeneity can be quantified with noninvasive imaging. We developed spatial habitats based on multiparametric physiologic MRI and evaluated associations between temporal changes in these habitats and progression-free survival (PFS) after concurrent chemoradiotherapy (CCRT) in patients with glioblastoma.

Experimental Design:

Ninety-seven patients with isocitrate dehydrogenase (IDH)-wildtype glioblastoma were enrolled and two serial MRI examinations after CCRT were analyzed. Cerebral blood volumes and apparent diffusion coefficients were grouped using k-means clustering into three spatial habitats. Associations between temporal changes in spatial habitats and PFS were investigated using Cox proportional hazard modeling. The performance of significant predictors for PFS and overall survival (OS) was measured using a discrete increase of habitat (habitat risk score) in a temporal validation set from a prospective registry (n = 53, ClinicalTrials.gov NCT02619890). The site of progression was matched with the spatiotemporal habitats.

Results:

Three spatial habitats of hypervascular cellular, hypovascular cellular, and nonviable tissue were identified. A short-term increase in the hypervascular cellular habitat (HR, 40.0; P = 0.001) and hypovascular cellular habitat was significantly associated with shorter PFS (HR, 3.78; P < 0.001) after CCRT. Combined with clinical predictors, the habitat risk score showed a C-index of 0.79 for PFS and 0.74 for OS and stratified patients with short, intermediate, and long PFS (P = 0.016). An increase in the hypovascular cellular habitat predicted tumor progression sites.

Conclusions:

Hypovascular cellular habitats derived from multiparametric physiologic MRIs may be useful predictors of clinical outcomes in patients with posttreatment glioblastoma.

Translational Relevance

Spatial heterogeneity in glioblastoma indicates poor clinical prognosis and drives resistance to therapy. The spatial habitats of posttreatment glioblastoma after concurrent chemoradiation therapy (CCRT) can be reflected by clustering of multiparametric physiologic MRI, including diffusion-weighted and perfusion-weighted MRI, which are widely available and can be repeated along the treatment course. Here, we define the three spatial habitats of hypervascular cellular tumor, hypovascular cellular tumor, and nonviable tissue within the contrast-enhancing lesion. Validation using a prospective cohort showed that an immediate increase in the hypovascular cellular tumor after CCRT was strongly associated with shorter progression-free survival. The hypovascular cellular tumor site matched the site of progression. Particularly, our spatiotemporal habitat analysis enables identification of early tumor progression in a more specific manner than the current response assessment for brain tumor and can potentially guide local therapies in patients with posttreatment glioblastoma.

Intratumoral heterogeneity in glioblastoma represents a biologically complex spatial variation in gene expression, histopathology, and macroscopic structure (1). Importantly, spatial heterogeneity frequently indicates a poor prognosis (2), because treatment response is not uniform across the tumor and resistance to therapy occurs in different tumor regions. The combination of genetic heterogeneity with metabolic differences are known to result in spatial habitats, which comprise subpopulations of cells preserving distinct characteristics (3–5). Tumor regions showing differences in metabolism, vascularity, and cellularity can be captured with physiologic MRI, including cerebral blood volume (CBV) and apparent diffusion coefficient (ADC) mapping, which provide data directly reflecting biologic information (6, 7).

Current tumor analysis using histograms (8) or radiomics analysis (9) focuses on quantifying the heterogeneity and complexity by calculating the relationship between voxels (5). However, the spatial distribution of voxels within the tumor are not considered, making biologic validation difficult (10). Identification of similar voxels differs in the underlying assumptions and methodology from that of radiomics, which may define multiple subregions at the voxel level with a common tumor biology (3, 5). This shows spatial correspondence with genomics and histopathology and provides clinically useful information for identifying subregions of tumor progression that may require biopsy, are resistant to current therapy, and may become a therapeutic target in the future.

Single-parametric changes, such as parametric response maps, have been used for predicting the treatment response of glioma (11–13). However, binary thresholds are not absolute and depend on the acquisition and analysis of images (5). In addition, tracking individual voxels is both extremely difficult and cannot directly identify tumor subregions based on multiparametric imaging parameters (5). A data-driven approach using voxel-wise clustering may identify functionally coherent tumor subregions based on multiparametric physiologic MRI (5, 14).

The pathophysiologic correlates of ADC and CBV are cell density and necrosis (15, 16), and vessel density (17, 18), respectively, but the measurement of a single quantitative parameter is limiting owing to tumor heterogeneity and treatment-related changes (19). We hypothesized that the temporal changes in spatial habits derived from multiparametric physiologic MRIs potentially provide information on both spatial and temporal heterogeneity in posttreatment glioblastoma. The technique does not require voxel-by-voxel matching for serial follow-ups and accounts for newly appeared or disappeared voxels by creating spatial habitats with clustering of similar ADC and CBV values. The purpose of this study was to develop multiparametric physiologic MRI-based spatiotemporal habitat analysis and to validate the association between temporal changes in these habitats with progression-free survival (PFS) after concurrent chemoradiotherapy (CCRT) in patients with glioblastoma.

Study population: training and validation sets

This retrospective clinical study was approved by the institutional review board (IRB) of Asan Medical Center (local approval number: 2020-0363) and was conducted in compliance with the U.S. Health Insurance Portability and Accountability Act regulations and the Declaration of Helsinki. For all the patients in the validation set, written and signed informed consent was obtained as per IRB-approval protocol. Figure 1 shows the inclusion process. A total of 150 patients with a histopathologic diagnosis of isocitrate dehydrogenase (IDH)-wildtype glioblastoma according to the World Health Organization classification 2016 (1), and a contrast-enhancing lesion (CEL) on MRI after CCRT were identified from the neuro-oncologic database of our institution between April 2016 and August 2018. Patients were enrolled according to the following inclusion criteria: (i) minimum age of 18 years; (ii) underwent standard treatment including surgery, CCRT, and adjuvant temozolomide; (iii) underwent conventional and advanced MRI, including diffusion-weighted imaging (DWI) and dynamic susceptibility contrast (DSC) imaging; (iv) had a measurable CEL of more than 1 × 1 cm2 on the first postchemoradiation examination; (v) clinicoradiologic follow-up with at least three successive imaging examinations or pathologic confirmation to monitor treatment response; and (vi) image quality with no artifacts affecting the analysis. Fifty-three patients did not meet the inclusion criteria and were excluded: 15 patients did not undergo DSC imaging, 3 had poor-quality imaging/perfusion data, and 35 did not undergo 1-year follow-up or more than three successive MRIs. Therefore, 97 consecutive patients [mean age, 58.1 years (range, 21–86); 48 (50.5%) women] were included in the study to determine the physiologic MRI-based spatial habitats.

An independent temporal validation set of IDH-wildtype glioblastoma was used to evaluate the true clinical performance. This independent dataset was from a prospective brain glioma registry (ClinicalTrials.gov identifier: NCT02619890) and included patients diagnosed and followed up between September 2018 and May 2020. Informed consent was obtained from all patients in a prospective brain glioma registry. The final validation set consisted of 53 patients [mean age, 57.8 ± 13.4 years (range, 32–74); 28 (52.8%) women].

MRI acquisitions

The training and validation set shared the same imaging acquisition protocol. The brain tumor imaging protocol was acquired on a 3-T scanner (Ingenia 3.0 CX, Philips Healthcare) and included conventional and advanced sequences, including T2-weighted imaging (T2WI), T2-weighted fluid-attenuated inversion recovery (FLAIR), and precontrast and postcontrast T1-weighted imaging (T1WI). The FLAIR and contrast-enhanced T1WI imaging parameters are shown in Supplementary Table S1.

The DWI parameters were as follows: repetition time (TR)/echo time (TE), 3,000/56 ms; diffusion gradient encoding, b = 0 and 1,000 seconds/mm2; field of view (FOV), 250 × 250 mm; matrix, 256 × 256; and slice thickness/gap, 5/2 mm. ADC images were calculated from b = 1,000 and b = 0 second/mm2 DWI images.

DSC imaging was acquired using a gradient-echo echo-planar imaging protocol. A preload of 0.01 mmol/kg gadoterate meglumine (Dotarem; Guerbet) was administered, followed by a dynamic bolus of a standard dose of 0.1 mmol/kg gadoterate meglumine delivered at a rate of 4 mL/second by an MRI-compatible power injector (Spectris; Medrad). The contrast material bolus was followed by injecting 20 mL saline at the same injection rate. The DSC imaging parameters included the following: TR/TE, 1,808/40 ms; flip angle, 35°; FOV, 24 × 24 cm; slice thickness/gap, 5/2 mm; matrix, 128 × 128; total acquisition time, 1 minute and 54 seconds. The dynamic acquisition was performed with a temporal resolution of 1.5 seconds, and 60 dynamics were acquired. DSC imaging was acquired with complete tumor volume coverage and the same section orientation as that used in conventional MRI. The MRI examinations performed at the first visit after CCRT (examination 1, Post-CCRT #1) and the second visit after CCRT (examination 2, Post-CCRT #2) were used in the analysis.

Mask segmentation and advanced image processing

For the three-dimensional contrast-enhanced T1WI (3D CE-T1WI) and FLAIR data, skull stripping was performed with an algorithm (https://github.com/MIC-DKFZ/HD-BET) optimized for processing heterogeneous MRI data with varying degrees of pathology or posttreatment alterations. For all patients, lesion segmentation masks were generated on the 3D CE-T1WI and FLAIR data using 3D UNet-based method (https://github.com/MIC-DKFZ/nnUNet; ref. 20) from the PyTorch package version 1.1 in Python 3.7 (www.python.org).

For the DSC analysis, pharmacokinetic map calculation was performed using Nordic ICE (NordicNeuroLab). Using the integrated DSC module, which incorporates a relative CBV (rCBV) leakage-correction algorithm and manual noise thresholding, the amount of blood in a given volume of tissue was obtained and expressed as mL per 100 mL tissue. The calculations using the Weisskoff–Boxerman method based on time-dependent deviation of the pixel-wise concentration–time curve from a reference curve were assumed to be unaffected by leakage (21). The rCBV maps were normalized according to the normal-appearing white matter, to create normalized CBV (nCBV) maps.

The 3D CE-T1WI images of each patient in the training dataset were coregistered and resampled into isometric-voxel sizes to evaluate changes across consecutive scans. Then, the FLAIR, nCBV, and ADC images were coregistered and resampled to the iso-voxel CE-T1WI images using rigid transformations with six degrees of freedom in SPM package (version 12, www.fil.ion.ucl.ac.uk/spm/). This step allowed continuous slices without gap and enabled tracking of habitats by voxel-wise analysis.

The final voxel classifications based on nCBV and ADC values were implemented using a k-means clustering module in the scikit-learn python package.

Multiparametric physiologic MRI-based spatiotemporal habitats

Optimization of number of clusters

All voxels from all segmented masks in the training set were first aggregated; then, all voxels in the ADC and nCBV maps were split into three clusters, which were functionally coherent subregions within a CEL. The individual voxels in each cluster were grouped on the basis of their similarities and differences using the k-means clustering algorithm with squared Euclidean distances between voxel intensities as the similarity metric. Because the optimal number of clusters in a dataset is an important issue in k-means clustering, 2, 3, 4, and 5 clusters were initially set. The aggregation of all voxels from all enhancing lesions in the primary cohort and the three clusters are demonstrated in Supplementary Fig. S1. The three clusters demonstrated both ADC and nCBV effects, while the fourth cluster was dominated by nCBV and the fifth cluster by ADC and displayed a narrow range of clusters. Finally, three clusters were chosen based on ADC and nCBV differences and to avoid overparameterized models (22).

Population-level clustering

Three clusters were set using two distinct feature maps: cluster 1 represented “hypervascular cellular tumor” with high CBV value and low ADC value; cluster 2 represented “hypovascular cellular tumor” with low CBV value and low ADC value; and cluster 3 represented “nonviable tissue” with low CBV value and high ADC value. All voxels were allocated into one of the three clusters and were displayed as spatial habitats in the original image space. The range for boundary of spatial habitats was 4.37–4.44 for nCBV and 150–187 (×10−6 mm2/s) for ADC.

Calculation of spatiotemporal habitats

The two consecutive MR examinations acquired at the first visit after CCRT (examination 1, Post-CCRT #1) and the second visit after CCRT (examination 2, Post-CCRT#2) were used in the analysis. First, changes in the numbers of voxels in the entire enhancing lesion and in each habitat were calculated between the two examinations (examination 2 − examination 1). Second, percentage of each habitat to the CEL volume was calculated (number of voxels in habitat 1/number of voxels in CEL volume) and change in the percentages between the two examinations (examination 2 − examination 1) was calculated.

Clinical predictors and outcome definition

Clinical, radiologic, and pathologic data were retrieved from medical records. Baseline characteristics including sex, age, Karnofsky performance score (KPS; binary, score >70 or ≤70), initial tumor volume, O6-methylguanine-DNA-methyltransferase (MGMT) promoter status, and extent of surgery (biopsy, partial resection, or gross total resection) were collected.

The primary endpoint was PFS. Tumor progression on MRI was assessed according to the Response assessment in neuro-oncology (RANO) criteria (23), using serial measures of the product of the two largest cross-sectional diameters. An objective response was defined as a complete or partial response on two consecutive MRI acquisitions with reduced or stable doses of corticosteroids. The consecutive clinicoradiologic diagnoses were made by consensus between a neuro-oncologist (J.H. Kim, with 26 years of experience in neuro-oncology practice) and a neuroradiologist (H.S. Kim, with 21 years of experience in neuro-oncologic imaging) after complete imaging and medical chart review (23). The overall survival (OS) was calculated from the day of diagnosis until the day of death, using the information obtained from the national health care data linked to our hospital.

Comparison with RANO assessment

RANO assessment was separately performed by two central reviewers (Y.-H. Kim and J.E. Park with 10 years of experience in neuro-oncologic practice and imaging, respectively), who reached consensus and had not participated in any other imaging review. The RANO assessment was based on postoperative imaging and examination 2 (Post-CCRT #2), and the C-index was compared with that from spatiotemporal habitat.

Analysis of site of progression

The site of progression was analyzed using the follow-up examination at the time of progression. The CEL volume at the time of progression was matched to the habitats at examination 2 (Post-CCRT #2). The overlap of each spatiotemporal habitat with the CEL volume at the time of progression was calculated using the DICE similarity coefficient, |{\rm{DICE}} = {\frac{{2| {P \cap R} |}}{{| P | + | R |}}$| (24). P indicates each spatiotemporal habitat and R indicates the CEL volume at the time of progression. The DICE ranges between 0 (no overlap) and 1 (perfect agreement).

Statistical analysis

Power calculation

The methodology written here adheres to the REMARK (Reporting Recommendations for Tumor Marker Prognostic Studies) recommendations (25). The sample size of the validation set was determined according to power analysis for a Cox regression of the log HR, which was performed using Power Analysis Software (PASS) version 15.0.7.

Spatiotemporal habitat to predict PFS

Univariate analysis was performed for analyzing the association of spatiotemporal habitats with PFS using Cox regression or the Kaplan–Meier method (log-rank test). HRs indicate relative change in hazard incurred by 1 unit increase in each parameter, and 20,000 voxels (20 k voxels) were considered as 1 unit in this study. For the significant predictors identified in the univariate Cox regression, an optimal cutoff for stratifying high- and low-risk groups was estimated using maxstat in R with 10-fold cross-validation, which ensured unbiased prediction within the sample (26). The z-score is the ratio of each regression coefficient to its SE, a Wald statistic which is asymptotically standard normal under the hypothesis that the corresponding β is zero (27). The z-score was used to demonstrate statistical significance of each spatiotemporal habitat.

Stratification of patients using spatiotemporal habitat

Using this cut-off from the exploratory analysis, a habitat risk score was developed to realize risk stratification for patients; each predictor was assigned a discrete score of 1 if it was larger than its cut-off, and 0, if it was lesser. For example, when two habitats relevant to PFS showed a discrete increase above the cut-off, the habitat risk score was 2.

Then, the C-index was calculated using independent predictors from multivariate Cox regression analyses to predict survival, by combining the habitat risk score and clinical predictors.

All results are reported as median with range or 95% confidence interval (CI) for continuous variables, and as frequencies or percentages for categorical variables. Differences in clinicopathologic factors between the training and validation cohorts were assessed using the χ2 test for discrete variables and Student t test for continuous variables. All statistical tests were two-sided, and a P value < 0.05 was considered statistically significant. Statistical analyses were performed using R software version 3.4.3 (https://www.r-project.org).

Demographic data of the study population

The clinical characteristics of the patients in the training and validation sets are summarized in Table 1. No differences in sex, age, initial KPS, treatment regimen, and MGMT promoter methylation status were noted between the training and validation cohorts.

Table 1.

Baseline clinical characteristics of the study patients.

Training setValidation set
Parameter(n = 97)(n = 53)P
Sex, n 
 Male/female 48/49 25/28 0.78 
Age, years 
 Mean ± SD (range) 58.2 ± 10.7 (21–86) 57.8 ± 13.4 (32–74) 0.06 
Initial tumor volume (mL) 
 Mean ± SD (range) 15.2 ± 9.6 (1.2–42.7) 14.7± 9.7 (2.6–31.9) 0.75 
Primary treatment, n (%)   0.25 
 Extent of resection    
  Gross–total resection 57 (58.8%) 28 (52.8%)  
  Subtotal resection 27 (27.8%) 21 (39.6%)  
  Biopsy 13 (13.4%) 4 (7.6%)  
 Adjuvant treatment   0.19 
  RT + TMZ 96 (98.9%) 52 (98.1%)  
  RT only 1 (1.1%) 1 (1.9%)  
KPS at treatment initiation, n (%)    
 >70 63 (64.9%) 40 (75.5%) 0.18 
 ≤70 34 (35.1%) 13 (24.5%)  
IDH mutation status, n (%)   >0.99 
 Negative 97 (100%) 53 (100%)  
MGMT promoter status, n (%)   0.48a 
 Methylated 21 (21.6%) 24 (45.3%)  
 Unmethylated 29 (29.9%) 25 (47.2%)  
 NA 47 (48.5%) 4 (7.5 %)  
Imaging interval between examination 1 and 2, days   0.052 
 Mean ± SD 75 ± 22 (24–116) 67 ± 22 (24–111)  
 Median follow-up time, years (range) 3.3 (1.1–4.5) 2.7 (1.4–3.6) 0.056 
Training setValidation set
Parameter(n = 97)(n = 53)P
Sex, n 
 Male/female 48/49 25/28 0.78 
Age, years 
 Mean ± SD (range) 58.2 ± 10.7 (21–86) 57.8 ± 13.4 (32–74) 0.06 
Initial tumor volume (mL) 
 Mean ± SD (range) 15.2 ± 9.6 (1.2–42.7) 14.7± 9.7 (2.6–31.9) 0.75 
Primary treatment, n (%)   0.25 
 Extent of resection    
  Gross–total resection 57 (58.8%) 28 (52.8%)  
  Subtotal resection 27 (27.8%) 21 (39.6%)  
  Biopsy 13 (13.4%) 4 (7.6%)  
 Adjuvant treatment   0.19 
  RT + TMZ 96 (98.9%) 52 (98.1%)  
  RT only 1 (1.1%) 1 (1.9%)  
KPS at treatment initiation, n (%)    
 >70 63 (64.9%) 40 (75.5%) 0.18 
 ≤70 34 (35.1%) 13 (24.5%)  
IDH mutation status, n (%)   >0.99 
 Negative 97 (100%) 53 (100%)  
MGMT promoter status, n (%)   0.48a 
 Methylated 21 (21.6%) 24 (45.3%)  
 Unmethylated 29 (29.9%) 25 (47.2%)  
 NA 47 (48.5%) 4 (7.5 %)  
Imaging interval between examination 1 and 2, days   0.052 
 Mean ± SD 75 ± 22 (24–116) 67 ± 22 (24–111)  
 Median follow-up time, years (range) 3.3 (1.1–4.5) 2.7 (1.4–3.6) 0.056 

Abbreviations: CCRT, concurrent chemoradiation therapy; IDH, isocitrate dehydrogenase; KPS, Karnofsky performance score; MGMT, O6-methylguanine-DNA-methyltransferase gene methylation status; NA, information not available; RT, radiotherapy; TMZ, temozolomide.

aχ2 test P value was calculated on the basis of available data.

The data from 42 patients (28%) who were alive at the time of analysis were censored and included. The median OS time of the deceased patients was 15 months (interquartile range, 15–29 months) in the training set and 13 months (interquartile range, 6–18 months) in the validation set. The median follow-up periods were 39 months and 19 months for patients in the training set and validation set, respectively.

The mean interval between the two serial examinations (Post-CCRT#1 and #2) was 75 ± 22 days in the training set (mean ± SD) and 67 ± 22 days in the validation set (mean ± SD); no significant difference was observed between the two datasets.

Spatiotemporal habitats associated with PFS in the training set

The results of the univariate analysis to evaluate the change in three spatiotemporal habitats are summarized in Table 2. The distribution of the average percentage of the three spatiotemporal habitats (hypervascular cellular cluster, hypovascular cellular cluster, and nonviable tissue cluster) was similar between examination 1 and 2; for examination 1, the values were 25.0%, 50.1%, and 25.9%, respectively, and for examination 2, 23.5%, 52.1%, and 24.4%, respectively, for each habitat (Supplementary Fig. S2). The habitat characteristics are shown in Supplementary Table S2.

Table 2.

Exploratory analysis of the spatial habitat to predict PFS in posttreatment IDH-wildtype glioblastoma.

Progression-free survival
Change between examinations 1 and 2 (examination 2 − examination 1)HR95% CIP (z statistics)
Spatiotemporal habitats (20 k voxel numbers) 
 Increase in hypervascular cellular habitat 40.00 4.71–339.2 0.001 
 Increase in hypovascular cellular habitat 3.78 2.38–6.01 <0.001 
 Increase in nonviable tissue habitat 2.59 0.29–23.2 0.39 
Spatiotemporal habitats (%) 
 Increase in percentage of hypervascular cellular habitat 1.01 0.97–1.04 0.72 
 Increase in percentage of hypovascular cellular habitat 1.01 1.00–1.02 0.02 
 Increase in nonviable tissue habitat 0.99 0.98–1.00 0.46 
Progression-free survival
Change between examinations 1 and 2 (examination 2 − examination 1)HR95% CIP (z statistics)
Spatiotemporal habitats (20 k voxel numbers) 
 Increase in hypervascular cellular habitat 40.00 4.71–339.2 0.001 
 Increase in hypovascular cellular habitat 3.78 2.38–6.01 <0.001 
 Increase in nonviable tissue habitat 2.59 0.29–23.2 0.39 
Spatiotemporal habitats (%) 
 Increase in percentage of hypervascular cellular habitat 1.01 0.97–1.04 0.72 
 Increase in percentage of hypovascular cellular habitat 1.01 1.00–1.02 0.02 
 Increase in nonviable tissue habitat 0.99 0.98–1.00 0.46 

Note: HRs reported here indicate the relative change in hazard that a 1-unit (20,000 voxels) increase in each imaging parameter incurs.

Among the various changes in spatial habitats, an increase in the hypovascular cellular habitat (20 k voxels) was most significantly associated with a relatively short PFS (HR, 3.78; 95% CI, 2.38–6.00; P < 0.001). An increase in the hypervascular cellular habitat (HR, 40.0; 95% CI, 4.71–339.2; P = 0.001; z-score = 3.38) and percentage of the hypovascular cellular habitat were also associated with a relatively short PFS (HR, 1.01; 95% CI, 1.00–1.02; P = 0.02; z-score = 2.31); however, the association was weaker than that observed for the increase in the hypovascular cellular habitat (z-score = 5.6). Figure 2 demonstrates the association of PFS with change in the hypovascular cellular habitat. An increase in the hypervascular cellular habitat or nonviable tissue habitat was not associated with PFS (P > 0.05).

Figure 1.

Flow diagram describing the process of patient inclusion.

Figure 1.

Flow diagram describing the process of patient inclusion.

Close modal
Figure 2.

A, Demonstration of the three spatial habitats defined by clustered voxels from normalized ADC and nCBV maps from a 55-year-old patient. The hypervascular cellular habitat (red color) shows high nCBV and low ADC, the hypovascular cellular habitat (green color) shows low nCBV and low ADC, and the nonviable tissue habitat (blue color) shows low nCBV and high ADC. B, Results of two serial examinations performed at 8 and 12 weeks after the completion of CCRT [examinations 1 (post-CCRT #1) and 2 (post-CCRT #2)] from a patient show an increase in the CEL volume and hypovascular cellular habitat, but a decrease in nonviable tissue habitat. The results of the confirmatory examination after 4 weeks show tumor progression.

Figure 2.

A, Demonstration of the three spatial habitats defined by clustered voxels from normalized ADC and nCBV maps from a 55-year-old patient. The hypervascular cellular habitat (red color) shows high nCBV and low ADC, the hypovascular cellular habitat (green color) shows low nCBV and low ADC, and the nonviable tissue habitat (blue color) shows low nCBV and high ADC. B, Results of two serial examinations performed at 8 and 12 weeks after the completion of CCRT [examinations 1 (post-CCRT #1) and 2 (post-CCRT #2)] from a patient show an increase in the CEL volume and hypovascular cellular habitat, but a decrease in nonviable tissue habitat. The results of the confirmatory examination after 4 weeks show tumor progression.

Close modal

Stratification based on change in the CEL volume and spatiotemporal habitat

The optimal cut-off of the spatiotemporal habitat needed for stratifying the longer and shorter PFS groups was an increase of more than 130 voxels in the hypovascular cellular habitat. The cut-off stratified the longer and shorter PFS groups with significant difference in the log-rank test (P = 0.016). Supplementary Figure S3 shows the Kaplan–Meier survival curve and the risk table generated on the basis of the hypovascular cellular habitat. The optimal cutoff for the hypervascular cellular habitat was an increase of more than 0 voxel. The habitat risk score, based on the discrete increase of the hypervascular and hypovascular cellular habitats was calculated for the patients. For example, when a patient showed an increase of both hypervascular cellular habitat (>0 voxel) and hypovascular cellular habitat (>130 voxels), the habitat risk score was 2.

In the training set, the habitat risk score stratified patients according to low risk (median PFS, 518 days; 95% CI, 441–672 days), intermediate risk (median PFS, 316 days; 95% CI, 180–527 days), and high risk (median PFS, 160 days; 95% CI, 149–279 days; log-rank test; P < 0.001).

A multivariate Cox analysis was conducted with the habitat risk score, age, KPS scores, extent of surgery, and the initial tumor volume as the factors. Only a high initial KPS score (HR, 0.43; 95% CI, 0.25–0.75; P = 0.003) and the habitat risk score (HR, 2.17; 95% CI, 1.59–2.96; P < 0.001) were identified as the independent predictors of PFS. The C-index values for PFS and OS were 0.79 (95% CI, 0.72–0.85) and 0.67 (95% CI, 0.64–0.71), respectively.

Independent validation in a temporal validation set

Power calculation for the validation set

The survival prediction value for our model was compared with that of the null model. The hypovascular cellular habitat showed HR 3.78, and HR 3.00 was conservatively set as the meaningful difference. The event rate (progression) was 0.7 and SD of the hypovascular cellular habitat was 0.9 in the training set. A sample size of 49 would achieve 80% power to have meaningful survival prediction based on HR 3.00, with an estimated SD of differences of 1.0 with a significance level (alpha) of 0.05 using a Cox regression model. Thus, the sample size of 53 in the validation set was sufficient for predicting survival.

The three spatial habitats displayed similar patterns in the validation cohort in contrast to those displayed in the training cohort. The habitat risk score successfully stratified patients according to low risk (median PFS, 329 days; 95% CI, 229–518 days), intermediate risk (median PFS, 122 days; 95% CI, 86–180 days), and high risk (median PFS, 80 days; 95% CI, 78–82 days; log-rank test; P = 0.02; Fig. 3). In combination with initial KPS score and the habitat risk score, the C-index values for PFS and OS were 0.79 (95% CI, 0.72–0.86) and 0.74 (95% CI, 0.70–0.78), respectively.

Figure 3.

Risk score estimated using the increase in both the hypervascular and hypovascular cellular habitats, based on which patients were stratified into high-, intermediate-, and low-risk groups in the validation cohort (log-rank test; P = 0.016).

Figure 3.

Risk score estimated using the increase in both the hypervascular and hypovascular cellular habitats, based on which patients were stratified into high-, intermediate-, and low-risk groups in the validation cohort (log-rank test; P = 0.016).

Close modal

Comparison with RANO-based assessment

The C-index for the RANO-based prediction was 0.58 (95% CI, 0.56–0.61) in the training set and 0.54 (95% CI, 0.50–0.58) in the validation set. The C-index from the habitat risk score only was 0.67 (95% CI, 0.64–0.71) in the training set and 0.71 (95% CI, 0.66–0.75) in the validation set, which was significantly better at predicting OS (P < 0.001) than that from RANO assessment.

Spatiotemporal habitats and site of progression

The tumor progressed in 46 patients (86.7%) in the validation set. Figure 4 demonstrates the relationship between the temporal changes in the hypovascular cellular habitat and the site of progression. The DICE index for the site of progression was calculated between examination 2 and the time of recurrence (examination 3). The mean DICE index was 0.13 (range, 0–0.6; SD 0.10) for hypervascular cellular habitat, 0.34 (range, 0–0.72; SD 0.12) for hypovascular cellular habitat, and 0.07 for nonviable tissue (range, 0–0.18; SD, 0.07). The DICE index for the hypovascular cellular habitat was highest in 34 of the 46 cases (73.9%).

Figure 4.

A representative example showing the relationship between the increase in the hypovascular cellular habitat and site of progression from a 62-year-old patient. A, The three spatial habitats in this patient defined by clustered voxels. B, Spatial mapping shows an increase in the hypovascular cellular habitat at the anterior portion of the enhancing masses, observed between 4 weeks (examination 1) and 16 weeks (examination 2) after the completion of CCRT. Subsequently, the hypovascular cellular habitat became the site of recurrence at the 3-month follow-up with a confirmatory exam.

Figure 4.

A representative example showing the relationship between the increase in the hypovascular cellular habitat and site of progression from a 62-year-old patient. A, The three spatial habitats in this patient defined by clustered voxels. B, Spatial mapping shows an increase in the hypovascular cellular habitat at the anterior portion of the enhancing masses, observed between 4 weeks (examination 1) and 16 weeks (examination 2) after the completion of CCRT. Subsequently, the hypovascular cellular habitat became the site of recurrence at the 3-month follow-up with a confirmatory exam.

Close modal

In this study, we demonstrated that PFS could be predicted in patients with IDH-wildtype glioblastoma using the spatiotemporal habitats, which were defined by multiparametric physiologic MRI values, quantified in the early stages of pathophysiologic alterations occurring after concurrent chemotherapy. An increase in the hypovascular cellular habitat with both low ADC and CBV values showed the most significant association with PFS. Furthermore, the increase in both the hypervascular and hypovascular cellular tumor was scored and stratified with a habitat risk score, and the patients in a prospective registry-based validation set were successfully stratified into groups with low, intermediate, and high risk of early progression.

The heterogeneous nature of glioblastomas after treatment (28) poses significant challenges in their evaluation using single quantitative parameters; such binary approaches are often affected by different image settings and algorithms (5). Furthermore, a threshold-based approach has an intrinsic limitation because it is difficult to label all tissue voxels when considerable heterogeneity is present, which results in tissue “gray-zones” (29). Alternatively, the clustering of voxels is advantageous because it requires no prior assumption (30) and is not affected by acquisitions that are normalized by the normal brain tissue. The centroids in k-means clustering can be stabilized without using a cut-off value (31); this approach has better implications for minimizing gray-zone tissue and applications in different imaging settings.

Temporal changes in physiologic MRI parameters when added to the information of spatial heterogeneity can enable clinical outcome prediction. Besides treatment-related changes such as low perfusion and high diffusion, the posttreatment glioblastoma tissues may also exhibit a spectrum of low-to-high perfusion and low diffusion (32). The hypervascular cellular habitat, with high CBV and low ADC values, represents the most malignant region of a tumor. However, the association between the change in the hypervascular cellular habitat with PFS was weaker than that of the hypovascular cellular habitat with PFS. Temporal changes in spatial habitats, such as “spatiotemporal habitats,” would be more helpful than a single snapshot of the posttreatment glioblastomas to predict PFS because of the inherent heterogeneity in the histopathology of posttreatment glioblastoma (33). Moreover, longitudinal analysis of an entire tumor with physiologic MRI might help to overcome this ambiguity and provide adequate information for subsequent treatment decisions.

The spatiotemporal habitat risk score was significantly better than the RANO-based assessment at predicting OS. This result is similar to that of a recent study based on the EORTC-26101 trial that showed the RANO assessment produced C-index of 0.57 (95% CI, 0.54–0.61) for predicting OS (34). Recent deep learning segmentation or T1-subtraction methods have emphasized the importance of objective tumor volume measurements for determining PFS (34, 35). As most glioblastomas recur locally (36, 37), a spatial habitat-driven approach may refine the volumetric analysis to predict their progression. From the results of our analysis, we speculate that ADC and CBV subdivides the contrast-enhancing tumor volume and this may provide more sensitive and specific information of tumor progression. Using the habitat risk score described herein, the change in spatiotemporal habitats can be easily determined, which might help to stratify patients with glioblastoma according to the risk of early progression.

Not only did the spatiotemporal habitats improve general survival performance, they also highlighted the importance of identifying subregions to determine intratumoral heterogeneity and identify regions with treatment resistance. We observed that the tumor recurrence foci matched the regions that showed short-term increase in the hypovascular cellular tumor habitat using quantitative analysis of the DICE index. The hypovascular cellular tumor habitat showed the highest DICE index with CEL at the time of progression. Parcellations ensure that the images are obtained at a single time point and that the functionally coherent regions can be followed up and matched to the tumor recurrence foci.

This study has several limitations. First, the study included a small cohort size, especially in the validation set. Second, strict pathologic correlations with image-based segmentations were lacking. However, such correlations would be very difficult to achieve, and unnecessary surgery should be avoided; therefore, we made assumptions regarding the clusters using a logical approach. Third, the three habitats defined here are relative indices of cellular habitats showing relatively low ADC values compared to those of other habitats. The b values selected in our study were only b = 0 and b = 1,000 seconds/mm2, thus, there is a risk of perfusion contamination in ADC values. One strategy is to obtain images with b value = 150–200 and 400–600 seconds/mm2 and use them to fit the data (38). This will provide more cellular components with less perfusion contamination, which is recommended in the future study. Fourth, although the study results were validated on a separate, prospectively enrolled cohort, this was a single-center study; thus, application of different protocols from other institutions with a relatively large number of patients is necessary. The normalization of quantitative values of ADC and CBV holds potential to increase reproducibility in a multicenter study. However, there is a big gap between taking the data as presented and turning them into a clinically useful tool that could guide patient therapy until it is proved in a multicenter setting, considering the cut-off for the habitat risk score depends on the training set. In the future, a prospective multicenter study using ADC and CBV in comparison with a RANO-based assessment to predict OS is further needed to confirm its utility. Finally, the technology used in this study is complex; the development of an easy-to-use software is warranted for widespread clinical adoption. A specialized postprocessing software for spatial habitat analysis and deep learning segmentation is being developed for clinical implementation.

In conclusion, the temporal change in spatial habitats derived from multiparametric physiologic MRIs may be a useful predictor of clinical outcome in patients with posttreatment glioblastoma. An increase in the hypovascular cellular habitat can be regarded as a strong imaging biomarker to predict early progression, as well the site of progression.

No disclosures were reported.

J.E. Park: Conceptualization, formal analysis, methodology, writing-original draft. H.S. Kim: Conceptualization, supervision, methodology, writing-review and editing. N. Kim: Software, formal analysis, investigation, visualization. S.Y. Park: Software, formal analysis. Y.-H. Kim: Data curation, project administration. J.H. Kim: Resources, project administration.

This research was supported by a National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIP; grant numbers: NRF-2020R1A2B5B01001707 and NRF-2020R1A2C4001748).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

1.
Louis
DN
,
Perry
A
,
Reifenberger
G
,
von Deimling
A
,
Figarella-Branger
D
,
Cavenee
WK
, et al
The 2016 World Health Organization classification of tumors of the central nervous system: a summary
.
Acta Neuropathol
2016
;
131
:
803
20
.
2.
Shipitsin
M
,
Campbell
LL
,
Argani
P
,
Weremowicz
S
,
Bloushtain-Qimron
N
,
Yao
J
, et al
Molecular definition of breast tumor heterogeneity
.
Cancer Cell
2007
;
11
:
259
73
.
3.
Dextraze
K
,
Saha
A
,
Kim
D
,
Narang
S
,
Lehrer
M
,
Rao
A
, et al
Spatial habitats from multiparametric MR imaging are associated with signaling pathway activities and survival in glioblastoma
.
Oncotarget
2017
;
8
:
112992
3001
.
4.
Lee
J
,
Narang
S
,
Martinez
J
,
Rao
G
,
Rao
A
. 
Spatial habitat features derived from multiparametric magnetic resonance imaging data are associated with molecular subtype and 12-month survival status in glioblastoma multiforme
.
PLoS One
2015
;
10
:
e0136557
.
5.
O'Connor
JPB
,
Rose
CJ
,
Waterton
JC
,
Carano
RAD
,
Parker
GJM
,
Jackson
A
. 
Imaging intratumor heterogeneity: role in therapy response, resistance, and clinical outcome
.
Clin Cancer Res
2015
;
21
:
249
57
.
6.
John
F
,
Bosnyak
E
,
Robinette
NL
,
Amit-Yousif
AJ
,
Barger
GR
,
Shah
KD
, et al
Multimodal imaging-defined subregions in newly diagnosed glioblastoma: impact on overall survival
.
Neuro Oncol
2019
;
21
:
264
73
.
7.
Kazerooni
AF
,
Nabil
M
,
Zadeh
MZ
,
Firouznia
K
,
Azmoudeh-Ardalan
F
,
Frangi
AF
, et al
Characterization of active and infiltrative tumorous subregions from normal tissue in brain gliomas using multiparametric MRI
.
J Magn Reson Imaging
2018
;
48
:
938
50
.
8.
Ryu
YJ
,
Choi
SH
,
Park
SJ
,
Yun
TJ
,
Kim
JH
,
Sohn
CH
. 
Glioma: application of whole-tumor texture analysis of diffusion-weighted imaging for the evaluation of tumor heterogeneity
.
PLoS One
2014
;
9
:
e108335
.
9.
Sala
E
,
Mema
E
,
Himoto
Y
,
Veeraraghavan
H
,
Brenton
JD
,
Snyder
A
, et al
Unravelling tumour heterogeneity using next-generation imaging: radiomics, radiogenomics, and habitat imaging
.
Clin Radiol
2017
;
72
:
3
10
.
10.
Park
JE
,
Park
SY
,
Kim
HJ
,
Kim
HS
. 
Reproducibility and generalizability in radiomics modeling: possible strategies in radiologic and statistical perspectives
.
Korean J Radiol
2019
;
20
:
1124
37
.
11.
Galban
CJ
,
Chenevert
TL
,
Meyer
CR
,
Tsien
C
,
Lawrence
TS
,
Hamstra
DA
, et al
Prospective analysis of parametric response map-derived MRI biomarkers: identification of early and distinct glioma response patterns not predicted by standard radiographic assessment
.
Clin Cancer Res
2011
;
17
:
4751
60
.
12.
Galban
CJ
,
Chenevert
TL
,
Meyer
CR
,
Tsien
C
,
Lawrence
TS
,
Hamstra
DA
, et al
The parametric response map is an imaging biomarker for early cancer treatment outcome
.
Nat Med
2009
;
15
:
572
6
.
13.
Tsien
C
,
Galban
CJ
,
Chenevert
TL
,
Johnson
TD
,
Hamstra
DA
,
Sundgren
PC
, et al
Parametric response map as an imaging biomarker to distinguish progression from pseudoprogression in high-grade glioma
.
J Clin Oncol
2010
;
28
:
2293
9
.
14.
Berens
ME
,
Sood
A
,
Barnholtz-Sloan
JS
,
Graf
JF
,
Cho
S
,
Kim
S
, et al
Multiscale, multimodal analysis of tumor heterogeneity in IDH1 mutant vs. wild-type diffuse gliomas
.
PLoS One
2019
;
14
:
e0219724
.
15.
Asao
C
,
Korogi
Y
,
Kitajima
M
,
Hirai
T
,
Baba
Y
,
Makino
K
, et al
Diffusion-weighted imaging of radiation-induced brain injury for differentiation from tumor recurrence
.
AJNR Am J Neuroradiol
2005
;
26
:
1455
60
.
16.
Lemasson
B
,
Galban
CJ
,
Boes
JL
,
Li
YH
,
Zhu
Y
,
Heist
KA
, et al
Diffusion-weighted MRI as a biomarker of tumor radiation treatment response heterogeneity: a comparative study of whole-volume histogram analysis versus voxel-based functional diffusion map analysis
.
Transl Oncol
2013
;
6
:
554
61
.
17.
Fatterpekar
GM
,
Galheigo
D
,
Narayana
A
,
Johnson
G
,
Knopp
E
. 
Treatment-related change versus tumor recurrence in high-grade gliomas: a diagnostic conundrum–use of dynamic susceptibility contrast-enhanced (DSC) perfusion MRI
.
AJR Am J Roentgenol
2012
;
198
:
19
26
.
18.
Barajas
RF
 Jr
,
Chang
JS
,
Segal
MR
,
Parsa
AT
,
McDermott
MW
,
Berger
MS
, et al
Differentiation of recurrent glioblastoma multiforme from radiation necrosis after external beam radiation therapy with dynamic susceptibility-weighted contrast-enhanced perfusion MR imaging
.
Radiology
2009
;
253
:
486
96
.
19.
Hygino da Cruz
LC
 Jr
,
Rodriguez
I
,
Domingues
RC
,
Gasparetto
EL
,
Sorensen
AG
. 
Pseudoprogression and pseudoresponse: imaging challenges in the assessment of posttreatment glioma
.
AJNR Am J Neuroradiol
2011
;
32
:
1978
85
.
20.
Isensee
F
,
Schell
M
,
Pflueger
I
,
Brugnara
G
,
Bonekamp
D
,
Neuberger
U
, et al
Automated brain extraction of multisequence MRI using artificial neural networks
.
Hum Brain Mapp
2019
;
40
:
4952
64
.
21.
Weisskoff
R
,
Boxerman
J
,
Sorensen
A
,
Kulke
S
,
Campbell
T
,
Rosen
B
.
Simultaneous blood volume and permeability mapping using a single Gd-based contrast injection. In: Proceedings of the Second Meeting of the Society of Magnetic Resonance
.
Berkeley, CA
:
International Society for Magnetic Resonance in Medicine
, 
1994
.
22.
Gull
SF
. 
Bayesian inductive inference and maximum entropy
.
In
:
Erickson
GJ
,
Smith
CR
,
editors
.
Maximum-entropy and bayesian methods in science and engineering: foundations
.
Dordrecht, Netherlands
:
Springer Netherlands
; 
1988
.
p.
53
74
.
23.
Wen
PY
,
Macdonald
DR
,
Reardon
DA
,
Cloughesy
TF
,
Sorensen
AG
,
Galanis
E
, et al
Updated response assessment criteria for high-grade gliomas: response assessment in neuro-oncology working group
.
J Clin Oncol
2010
;
28
:
1963
72
.
24.
Dice
LR
. 
Measures of the amount of ecologic association between species
.
Ecology
1945
;
26
:
297
302
.
25.
McShane
LM
,
Altman
DG
,
Sauerbrei
W
,
Taube
SE
,
Gion
M
,
Clark
GM
, et al
REporting recommendations for tumour MARKer prognostic studies (REMARK)
.
Br J Cancer
2005
;
93
:
387
91
.
26.
Ingrisch
M
,
Schneider
MJ
,
Norenberg
D
,
Negrao de Figueiredo
G
,
Maier-Hein
K
,
Suchorska
B
, et al
Radiomic analysis reveals prognostic information in T1-weighted baseline magnetic resonance imaging in patients with glioblastoma
.
Invest Radiol
2017
;
52
:
360
6
.
27.
Fox,
J
. 
Cox proportional-hazards regression for survival data. An R and S-PLUS companion to applied regression; 2002
. Available from: https://socialsciences.mcmaster.ca/jfox/Books/Companion-2E/appendix/Appendix-Cox-Regression.pdf.
28.
Jung
V
,
Romeike
BF
,
Henn
W
,
Feiden
W
,
Moringlane
JR
,
Zang
KD
, et al
Evidence of focal genetic microheterogeneity in glioblastoma multiforme by area-specific CGH on microdissected tumor cells
.
J Neuropathol Exp Neurol
1999
;
58
:
993
9
.
29.
Park
JE
,
Kim
HS
,
Goh
MJ
,
Kim
SJ
,
Kim
JH
. 
Pseudoprogression in patients with glioblastoma: assessment by using volume-weighted voxel-based multiparametric clustering of MR imaging data in an independent test set
.
Radiology
2015
;
275
:
792
802
.
30.
Kanungo
T
,
Mount
DM
,
Netanyahu
NS
,
Piatko
CD
,
Silverman
R
,
Wu
AY
. 
An efficient k-means clustering algorithm: analysis and implementation
.
IEEE Trans Pattern Anal
2002
;
24
:
881
92
.
31.
Wu
M-N
,
Lin
C-C
,
Chang
C-C
. 
Brain tumor detection using color-based k-means clustering segmentation
.
IEEE
2007
;
2
:
245
50
.
32.
Prager
AJ
,
Martinez
N
,
Beal
K
,
Omuro
A
,
Zhang
Z
,
Young
RJ
. 
Diffusion and perfusion MRI to differentiate treatment-related changes including pseudoprogression from recurrent tumors in high-grade gliomas with histopathologic evidence
.
AJNR Am J Neuroradiol
2015
;
36
:
877
85
.
33.
Melguizo-Gavilanes
I
,
Bruner
JM
,
Guha-Thakurta
N
,
Hess
KR
,
Puduvalli
VK
. 
Characterization of pseudoprogression in patients with glioblastoma: is histology the gold standard?
J Neurooncol
2015
;
123
:
141
50
.
34.
Kickingereder
P
,
Isensee
F
,
Tursunova
I
,
Petersen
J
,
Neuberger
U
,
Bonekamp
D
, et al
Automated quantitative tumour response assessment of MRI in neuro-oncology with artificial neural networks: a multicentre, retrospective study
.
Lancet Oncol
2019
;
20
:
728
40
.
35.
Ellingson
BM
,
Cloughesy
TF
,
Lai
A
,
Nghiemphu
PL
,
Mischel
PS
,
Pope
WB
. 
Quantitative volumetric analysis of conventional MRI response in recurrent glioblastoma treated with bevacizumab
.
Neuro Oncol
2011
;
13
:
401
9
.
36.
Sherriff
J
,
Tamangani
J
,
Senthil
L
,
Cruickshank
G
,
Spooner
D
,
Jones
B
, et al
Patterns of relapse in glioblastoma multiforme following concomitant chemoradiotherapy with temozolomide
.
Br J Radiol
2013
;
86
:
20120414
.
37.
Adeberg
S
,
Konig
L
,
Bostel
T
,
Harrabi
S
,
Welzel
T
,
Debus
J
, et al
Glioblastoma recurrence patterns after radiation therapy with regard to the subventricular zone
.
Int J Radiat Oncol Biol Phys
2014
;
90
:
886
93
.
38.
Le Bihan
D
,
Breton
E
,
Lallemand
D
,
Aubin
ML
,
Vignaud
J
,
Laval-Jeantet
M
. 
Separation of diffusion and perfusion in intravoxel incoherent motion MR imaging
.
Radiology
1988
;
168
:
497
505
.