Prospective longitudinal analysis of imaging-based spatiotemporal tumor habitats in glioblastoma, IDH-wild type: implication in patient outcome using multiparametric physiologic MRI

Abstract

Background

Physiologic MRI-based tumor habitat analysis has the potential to predict patient outcomes by identifying the spatiotemporal habitats of glioblastoma. This study aims to prospectively validate the cut-off for tumor progression obtained from tumor habitat analysis based on physiologic MRI in ascertaining time-to-progression (TTP) and the site of progression in glioblastoma patients following concurrent chemoradiotherapy (CCRT).

Methods

In this prospective study (ClinicalTrials.gov ID: NCT02613988), we will recruit patients with IDH-wild type glioblastoma who underwent CCRT and obtained immediate post-operative and three serial post-CCRT MRI scans within a three-month interval, conducted using diffusion-weighted imaging and dynamic susceptibility contrast imaging. Voxels from cerebral blood volume and apparent diffusion coefficient maps will be grouped using k-means clustering into three spatial habitats (hypervascular cellular, hypovascular cellular, and nonviable tissue). The spatiotemporal habitats of the tumor will be evaluated by comparing changes in each habitat between the serial MRI scans (post-operative and post-CCRT #1, #2, and #3). Associations between spatiotemporal habitats and TTP will be analyzed using cox proportional hazard modeling. The site of progression will be matched with spatiotemporal habitats.

Discussion

The perfusion- and diffusion-derived tumor habitat in glioblastoma is expected to stratify TTP and may serve as an early predictor for tumor progression in patients with IDH wild-type glioblastoma.

Trial registration

ClinicalTrials.gov ID: NCT02613988.

Background

Glioblastoma, known for its intratumoral heterogeneity, exhibits complex spatial variations in gene expression, histopathology, and macroscopic structure [1]. This heterogeneity is associated with a poor prognosis due to varied treatment responses and the development of therapy resistance in different tumor regions [2]. Physiological MRI techniques, such as cerebral blood volume (CBV) and apparent diffusion coefficient (ADC) mapping, enable the detection of distinct tumor regions that exhibit variations in metabolism, vascularity, and cellularity [3]. ADC provides insights into cell density and necrosis [4], while CBV correlates with vessel density [5]. While several imaging techniques, such as histograms, texture analysis, and radiomics, have been employed to quantify intratumoral heterogeneity based on imaging parameters, these methods often overlook spatial information, preventing the grouping of similar voxels [6, 7].

To overcome this limitation, tumor habitat analysis has emerged as a promising approach to distinguish subregions within a heterogeneous tumor by identifying voxels that share common tumor biology [8]. By utilizing voxel-wise clustering, this novel approach identifies multiple subregions within a tumor that share common tumor biology, offering valuable clinical insights into tumor subregions associated with progression, therapy resistance, and potential therapeutic targets [9]. Thus, clustering of multiparametric physiologic MRI, including diffusion-weighted and perfusion-weighted MRI, can reflect spatial habitats of post-treatment glioblastoma following concurrent chemoradiation therapy (CCRT). Furthermore, analyzing the temporal changes in spatial habitats derived from multiparametric physiologic MRIs could provide valuable insights into both spatial and temporal heterogeneity in post-treatment glioblastoma.

Recent research [10] has identified three spatial habitats within post-treatment glioblastoma: hypervascular cellular tumor, hypovascular cellular tumor, and nonviable tissue. Among these spatial habitats, an immediate increase in the hypovascular cellular tumor domain was strongly associated with poor prognosis following CCRT. Further, the localization of the hypovascular cellular tumor site also correlated with the site of disease progression. However, the true clinical performance of physiological MRI-based tumor habitat analysis in predicting patient outcomes has not been demonstrated in prospective studies. It is also necessary to compare this approach with the widely adopted RANO-based assessment for predicting time-to-progression (TTP).

Our ongoing clinical study therefore aims to prospectively validate physiologic MRI-based tumor habitat analysis in predicting TTP and identifying the site of progression after CCRT in patients with IDH-wild type glioblastoma.

Methods

Study design

This study is designed as a single-center, open-label, single-arm trial conducted at Asan Medical Center, a university-affiliated 2700-bed tertiary medical center located in Seoul, Republic of Korea. Enrollment of participants will continue until December 2025, with eligible participants being recruited from the hospital. The study protocol has received approval from the Institutional Review Board (IRB) of Asan Medical Center (2019 − 1259), and written informed consent will be obtained from each participant prior to enrollment. The imaging protocol was previously released as Early Response Assessment Using on 3T Advanced MR Imaging as Predictor of Long-term Treatment Response in Newly Diagnosed Glioblastomas (NCT02613988). This will be a validation study for imaging biomarker (tumor habitat analysis) obtained from a prospective imaging protocol enrollment from NCT02613988, utilizing voxel-wise clustering method of DWI and DSC imaging to represent tumor cellularity and vascularity with high generalizability.

Eligibility criteria

All study participants will be required to meet the following inclusion criteria and provide informed consent. Patients who do not meet the inclusion criteria will be excluded from the study. The flow of participant inclusion is schematically shown in Fig. 1.

Fig. 1

Flow chart of the study. Abbreviations: CCRT, chemoradiotherapy; TMZ, temozolomide

Inclusion criteria

(1) Adult patients aged over 18 years.

(2) Patients who were histologically diagnosed with IDH wild-type glioblastoma according to the World Health Organization classification 2021 [11].

(3) Patients scheduled for current standard treatment: maximum safe surgical resection followed by concurrent TMZ (75 mg/m²/day for 6 weeks) and RT (60 Gy in 30 fractions) and then six maintenance cycles of TMZ (150–200 mg/m²/day for the first 5 days of a 28-day cycle—TMZ) [12].

(4) Patients who provide written informed consent prior to any study-specific procedures.

Exclusion criteria

(1) Patients unable to receive MRI brain scan.

(2) Patients with no evidence of a measurable contrast-enhancing lesion of more than 1 × 1 cm² on the first post-CCRT examination.

MRI acquisitions

The brain tumor imaging protocol will utilize a 3-T scanner (Ingenia 3.0 CX, Philips Healthcare) and will include both structural and physiologic sequences: T2-weighted imaging (T2WI), fluid-attenuated inversion recovery (FLAIR) imaging, T1-weighted imaging (T1WI), DWI, DSC perfusion imaging, and contrast-enhanced (CE) T1WI.

The DWI parameters are as follows: repetition time (TR)/echo time (TE) at 3,000/56 ms, diffusion gradient encoding with b values of 0 and 1,000 s/mm², field of view (FOV) of 250 × 250 mm, matrix dimensions of 256 × 256, and a slice thickness/gap of 5/2 mm. ADC images will be calculated based on the DWI images acquired with b values of 1,000 and 0 s/mm².

DSC imaging will be conducted using a gradient-echo echo-planar imaging protocol. A preload of 0.01 mmol/kg gadoterate meglumine (Dotarem; Guerbet) will be administered, followed by a dynamic bolus of a standard dose of 0.1 mmol/kg gadoterate meglumine at a rate of 4 mL/second using an MRI-compatible power injector (Spectris; Medrad). Subsequently, 20 mL of saline will be injected at the same rate. The DSC imaging parameters are set as follows: TR/TE of 1,808/40 ms, a flip angle of 35°, FOV of 24 × 24 cm, slice thickness/gap of 5/2 mm, matrix dimensions of 128 × 128, and a total acquisition time of 1 min and 54 s. The dynamic acquisition will have a temporal resolution of 1.5 s, capturing a total of 60 dynamics. The DSC imaging will cover the entire tumor volume with the same section orientation as conventional MRI.

Reference standard for final diagnosis and endpoints

A final diagnosis of pseudoprogression (treatment-related change) or true progression will be confirmed pathologically by surgery when clinically indicated. In cases where second-look operations are not performed, consecutive clinicoradiological diagnoses will be established by consensus between two experts: 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. These diagnoses will be based on the RANO criteria [13]. Progression will be defined as the occurrence of any new lesion outside the radiation field or a gradual increase in the size of the contrast-enhancing lesion observed in more than two subsequent follow-up MRI examinations conducted at 2–3-month intervals, necessitating a prompt change in treatment. A diagnosis of progression will not be made if the increase was caused by comorbid conditions or concurrent medication). For enlarged contrast-enhancing lesions, the size criterion was an increase of more than 25% in the size of a measurable (> 1 cm) enhancing lesion according to the sum of the products of the perpendicular dimensions between the pre-CCRT and first post-CCRT MR images. This was modified from the RANO criteria for progressive disease [13] less than 12 weeks after CCRT completion.

On the other hand, pseudoprogression will be defined when no change in treatment is required at least 6 months after the end of CCRT. This diagnosis allows for a mild increase in contrast-enhancing lesions, as long as there is no treatment change during this time period. The cases of pseudoprogression will be classified as “non-progression”. The time of first progression will be calculated as the date of true progression for these cases.

The primary endpoint of the study is TTP, calculated from the day of initial diagnosis to the day of first documented progression. The secondary endpoint is the site of progression.

Mask segmentation and advanced image processing

To process the three-dimensional contrast-enhanced T1-weighted imaging (3D CE T1WI) and FLAIR data, a skull stripping algorithm optimized for heterogeneous MRI data will be employed (https://github.com/MIC-DKFZ/HD-BET). Lesion segmentation masks will be generated using a 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 will be performed using Nordic ICE (NordicNeuroLab). The integrated DSC module incorporates a relative CBV (rCBV) leakage correction algorithm and manual noise thresholding to quantify the amount of blood in a given volume of tissue, expressed as mL per 100 mL tissue. The Weisskoff-Boxerman method, which calculates pixel-wise concentration-time curve deviations from a reference curve, will be used for calculations assumed to be unaffected by leakage [14]. The rCBV maps will be normalized based on the normal-appearing white matter to create normalized CBV (nCBV) maps.

To analyze changes across consecutive scans, the 3D CE T1WI images obtained from each patient will be coregistered and resampled to have isometric-voxel sizes. Subsequently, FLAIR, nCBV, and ADC images will be coregistered and resampled to the iso-voxel CE T1WI images using rigid transformations with six degrees of freedom in the SPM package (version 12, www.fil.ion.ucl.ac.uk/spm/). This step ensures continuous slices without gaps and facilitates voxel-wise analysis for tracking habitats. The final voxel classifications based on nCBV and ADC values will be implemented using a k-means clustering module in the scikit-learn python package.

Multiparametric physiologic MRI-based spatiotemporal habitat analysis

Population-level clustering based on previous research

Using two distinct feature maps, three clusters will be established: cluster 1 representing the “hypervascular cellular tumor” characterized by a high CBV value and low ADC value, cluster 2 representing the “hypovascular cellular tumor” with low CBV value and low ADC value, and cluster 3 representing the “nonviable tissue” exhibiting low CBV value and high ADC value. The range for the boundary of the pre-trained and retrospectively validated spatial physiologic habitats was previously reported as 4.37–4.44 for nCBV and 150–187 (×10^− 6 mm²/s) for ADC [10].

Calculation of spatiotemporal habitats

The analysis will utilize four consecutive MR examinations: immediate post-operation (examination 1, post-op), the first visit after CCRT (examination 2, post-CCRT #1), the second visit after CCRT (examination 3, post-CCRT #2), and the third visit after CCRT (examination 4, post-CCRT #3), with a 3-month interval between each examination. Firstly, changes in the number of voxels within the entire enhancing lesion and within each habitat will be calculated between sequential examinations. Secondly, the percentage of each habitat relative to the contrast-enhancing lesion (CEL) volume will be calculated as the number of voxels in habitat 1 divided by the number of voxels in the CEL volume. We will then calculate changes in these percentages between sequential examinations.

Analysis of site of progression

The site of progression will be analyzed based on the follow-up examination at the time of progression. The volume of CEL at the time of progression will be matched with the habitats identified in examination 2 (post-CCRT#1) and 3 (post-CCRT #2). The overlap between each spatiotemporal habitat and the CEL volume at the time of progression will be quantified using the DICE similarity coefficient, DICE = $\:\frac{2\left|\text{P}\cap\:\text{R}\right|}{\left|\text{P}\right|+\left|\text{R}\right|}$ [15]. P represents each spatiotemporal habitat and R represents the CEL volume at the time of progression. The DICE ranges between 0 (no overlap) and 1 (perfect agreement).

The entire overview of our study process is depicted in Fig. 2. Representative cases of tumor habitat analysis is shown in Fig. 3.

Fig. 2

Process of spatiotemporal tumor habitat analysis using K-means clustering for ADC and CBV, and prediction of time-to-progression and site of progression. Abbreviations: CEL, contrast enhancing lesion; ADC, apparent diffusion coefficient; CBV, cerebral blood volume; CCRT, chemoradiotherapy

Fig. 3

Representative case of tumor habitat analysis from a 64-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. (A) Spatial mapping shows an increase in both hypervascular and hypovascular cellular habitat in post-CCRT examinations. The confirmatory scan after 4 weeks indicates tumor progression. (B) Spatial habitats defined by clustered voxels. Between postop and post-CCRT #1 examinations, the hypervascular cellular habitat increased by 104 voxels, and the hypovascular cellular habitat by 166 voxels, resulting in a habitat risk score of 2 points. Between post-CCRT first and second scans, hypervascular cellular habitat increased by 332 voxels, hypovascular cellular habitat by 1256 voxels, resulting in a habitat risk score of 2 points. Abbreviations: nCBV, normalized cerebral blood volume; ADC, apparent diffusion coefficient; CCRT, concurrent chemoradiotherapy

Statistical analysis

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 will be analyzed using descriptive statistics.

Sample sizecalculation

The methodology written here adheres to the REMARK (Reporting Recommendations for Tumor Marker Prognostic Studies) recommendations [16]. The sample size was calculated according to the Cox proportional hazards regression model with nonbinary covariates [17]. An event rate (progression) of 0.5 ~ 1.0 represent 51 progressions, with a sample size between 51 and 102 patients yielding expected power of 80% and an alpha error of 5%. We aim to recruit 100 patients over a 3-year period, with a 2-year follow-up for clinical assessments.

Use of spatiotemporal habitat to predict TTP

To analyze the association between spatiotemporal habitats and TTP, univariable analysis will be conducted using Cox proportional hazard regression or the Kaplan-Meier method (log-rank test). Hazard ratios indicate relative change in hazard incurred by 1 unit increase in each parameter; 20,000 voxels (20k voxels) defined a single unit in a previous study [10]. Additionally, a multivariable Cox analysis will be performed, considering the spatiotemporal habitat, age, KPS scores, extent of surgery, and initial tumor volume as factors predictive of TTP.

Risk stratification using spatiotemporal habitat

The habitat risk score, based on the discrete increase of the hypervascular and hypovascular cellular habitats will be calculated for the patients. The cutoff of a discrete score was calculated previously [10]: when a patient showed an increase of both hypervascular cellular habitat (> 0 voxel) and hypovascular cellular habitat (> 130 voxel), the habitat risk score was 2. Habitat risk score stratification between high-, intermediate-, and low-risk patients will be calculated using log-rank test. The log-rank test with p-value less than 0.05 would demonstrate that the risk score obtained from tumor habitat analysis successfully stratified patients with early progression (high-risk group).

All statistical analyses will be conducted using the R Statistical Package (version 3.6.3, Institute for Statistics and Mathematics, http://www.R-project.org). A p-value less than 0.05 will be considered statistically significant.

Discussion

This study aims to prospectively validate the predictive utility of physiological MRI-based tumor habitats in determining TTP and identifying the site of progression in patients with IDH-wild type glioblastoma after CCRT.

The inherent heterogeneity of post-treatment glioblastomas makes assessment using solitary quantitative parameters challenging [18]. Apart from changes induced by treatment, such as decreased perfusion and increased diffusion, the post-treatment glioblastoma tissues may demonstrate a spectrum of perfusion levels from low to high, along with decreased diffusion [19]. In this study, we hypothesize that ADC and CBV can more distinctly define the characteristics of glioblastoma after treatment, providing sensitive and specific information on tumor progression. Additionally, this approach may facilitate quantification of the initial pathophysiologic changes in post-treatment glioblastoma. Our previous research [10] categorized tumors into three groups: hypervascular cellular tumor, hypovascular cellular tumor, and nonviable tumor. An increase in the hypovascular cellular habitat correlated with a poor patient outcome, indicating its potential as a reliable imaging biomarker for early tumor progression prediction.

In the context of recurrent glioblastoma, re-resection is limited [20], and the presence of spatial heterogeneity in post-treatment glioblastoma hinders the achievement of adequate lesion sampling for effective histological analysis. Additionally, even with sufficient tissue sampling, explicit standards for histologic diagnosis of pseudoprogression, residual glioma, and recurrent glioma are currently lacking [21]. Therefore, we believe that utilizing advanced brain tumor imaging for radiographic guidance could facilitate optimal tissue sampling in cases of recurrent glioblastoma. Furthermore, if tumor habitat analysis enables non-invasive identification of the site of progression, its utility is anticipated to be significant.

In this study, we will analyze immediate post-operative and post-CCRT physiologic MRI scans in patients with IDH-wild type glioblastoma. Voxels from CBV and ADC maps will be grouped using k-means clustering into three spatial habitats. We will examine the temporal changes in these clusters and analyze the correlation between spatiotemporal habitats and patient outcomes.

In conclusion, we believe that tumor habitat analysis has the power to capture details of the spatiotemporal habitats of tumor vascularity and cellularity. This enhanced insight into the tumor environment has the potential to stratify TTP and serve as a useful predictor for early tumor progression and clinical outcomes in post-treatment glioblastoma patients.