Low perfusion compartments in glioblastoma quantified by advanced magnetic resonance imaging and correlated with patient survival

Highlights • Multiparametric MRI may identify two low perfusion compartments in glioblastoma.• The volume relates to better outcome, while the lactate relates to worse outcome.• Low perfusion compartment with restricted diffusivity may relates to resistance.

Glioblastoma is a highly aggressive primary brain malignancy in adults [1]. Despite advances in treatment, the median overall survival (OS) of patients remains low at 14.6 months [2]. Inconsistent response is a major challenge in treatment and could be caused by the extensive heterogeneity of this malignancy. Many genetically distinct cell populations can exist in the same tumor and display diverse treatment response [3,4].
One of the most fundamental traits of glioblastoma is the tumor-related angiogenesis and elevated perfusion, associated with a more invasive phenotype [5]. However, a potent angiogenesis inhibitor failed to demonstrate consistent benefits in clinical trials of de novo glioblastoma [6]. One possible explanation is the profound intratumor perfusion heterogeneity in glioblastomas, which is due to aberrant microvasculature and inefficient nutrient delivery. This heterogeneity can give rise to regions within tumors where the demand and supply of nutrients is mismatched [7]. Consequently, the sufficiently perfused habitats may hold the advantages for progression and proliferation, whereas the insufficiently perfused habitats may harbor a more acidic microenvironment than other tumor sub-regions [8], which may preferentially induce resistant clones to adjuvant therapy [9]. There is a rising need to understand the function of low perfusion compartments and evaluate their effects on treatment resistance.
Current clinical practice infers the low perfusion regions as the non-enhancing regions within contrast enhancement on postcontrast images, which can lead to non-specific results using conventional weighted images [10,11]. Recent studies suggested that quantitative imaging features are useful in reflecting intratumor habitats and tumor microenvironment [12,13]. As such, multiparametric imaging may allow a more comprehensive evaluation of the low perfusion compartments compared to the morphological heterogeneity visualized by structural magnetic resonance imaging (MRI).
The purpose of this study was to propose a method for quantifying low perfusion compartments in glioblastoma using multiparametric MRI and habitat imaging, and investigate their effects on treatment response and patient outcome. Our hypothesis is that multiparametric MRI may facilitate the identification of clinically relevant intratumor habitats that correlate with patient prognosis.
We integrated perfusion, diffusion and MR spectroscopic imaging with conventional imaging in our study. The relative cerebral blood volume (rCBV) calculated from perfusion imaging measures tumor vascularity [14]. The apparent diffusion coefficient (ADC) calculated from diffusion imaging provides information about tumor microstructure by measuring diffusivity of water molecules [15]. We used a thresholding method to visualize two low perfusion compartments with high/low diffusivity. Thus, two low perfusion compartments, distinguished by two potentially distinct properties, were visualized: one compartment with restricted diffusivity that may represent a sub-region with more microstructure adapting to the acidic conditions [12], and one compartment with increased diffusivity that may represent a sub-region with diminishing microstructure. We studied the metabolic signatures in each compartment using MR spectroscopy. Using multivariate survival analysis, we demonstrated that the volume and lactate level of these two compartments are clinically important.

Patient cohort
Patients with a radiological diagnosis of primary supratentorial glioblastoma suitable for maximal safe surgical resection were prospectively recruited from July 2010 to April 2015. All patients had a good performance status (World Health Organization performance status 0-1). Exclusion criteria included history of previous cranial surgery or radiotherapy/chemotherapy, or inability to undergo MRI scanning. This study was approved by the local institutional review board. Signed informed consent was obtained from all patients.
A total of 131 patients were recruited for the imaging scanning. After surgery, non-glioblastomas were excluded and 112 patients were included for regions of interest (ROI) analysis. Subgroups of patients with available MRS data (58 patients), DTI invasive phenotype data (64 patients) and survival data (80 patients) were analyzed. Patient recruitment and subgroups were summarized in Supplementary Fig. 1.

Treatment and response evaluation
All patients were on stable doses (8 mg/day) of dexamethasone. Tumor resection was performed with the guidance of neuronavigation (StealthStation, Medtronic) and 5-aminolevulinic acid fluorescence with other adjuvants, e.g., cortical and subcortical mapping, awake surgery, and intraoperative electrophysiology, to allow for maximal safe resection. Extent of resection was assessed according to the postoperative MRI scans within 72 hours as complete resection, partial resection of contrast-enhancing tumor or biopsy [16].
Patients received adjuvant therapy post-operatively according to their performance status. All patients were followed up according to the criteria of response assessment in neuro-oncology (RANO) [17], incorporating clinical and radiological criteria. Survivals were analyzed retrospectively in some cases when pseudoprogression was suspected if new contrast enhancement appeared within first 12 weeks after completing chemoradiotherapy.

Image processing
For each subject, all images were co-registered to T2-weighted images with an affine transformation, using the linear image registration tool functions [18] in Oxford Centre for Functional MRI of the Brain Software Library (FSL) v5.0.0 (Oxford, UK) [19].
DSC was processed and leakage correction was performed using NordicICE (NordicNeuroLab, Bergen, Norway). The arterial input function was automatically defined and rCBV was calculated. DTI images were processed with the diffusion toolbox in FSL [20], during which normalization and eddy current correction were performed. The isotropic component (p) and anisotropic component (q) were calculated using previous equations [21].
2D 1 H-MRSI were processed using LCModel (Provencher, Oakville, Ontario) and the concentrations of lactate (Lac) and macromolecule and lipid levels at 0.9 ppm (ML9) were calculated as a ratio to creatine (Cr). The quality and reliability of MRSI, and all spectra were evaluated using previous criteria [22].

Regions of interest and volumetric analysis
Tumor Regions of interest (ROIs) were manually segmented in 3D slicer v4.6.2 (https://www.slicer.org/) [23] by a neurosurgeon (ÂÂ, >7 years of experience) and a researcher (ÂÂ, >4 years of brain tumor image analysis experience) and reviewed by a neuroradiologist (ÂÂ, >8 years of experience) on post-contrast T1W and FLAIR images using a consistent criteria in each patient, and then cross-validated by comparing the similarity of the delineation using Dice similarity coefficient scores, blinded to the patient outcomes. The contrast-enhancing (CE) ROI was defined as all abnormalities within the contrast-enhancing rim on the post-contrast T1W images. FLAIR ROI was defined as the hyperintense abnormalities on FLAIR images. The interrater variability was tested using Dice similarity coefficient scores. For each individual patient, ROIs of normal-appearing white matter (NAWM) were manually segmented from the contralateral white matter as normal control. This region was typically 10 mm in diameter and located in the white matter which was furthest in distance from the tumor location and has no perceivable abnormalities [24]. All parametric maps of ADC and rCBV were normalized by the mean value in NAWM.
ADC-rCBV ROIs were generated using quartile values in MATLAB (v2016a, The MathWorks, Inc., Natick MA) from the CE ROI. The procedure is illustrated in Fig. 1. Firstly, ADC and rCBV values were obtained from each voxel within the CE ROI and pooled together as described previously [24]. The lowest quartile of the pooled rCBV values (rCBV L ) were interpreted as low perfusion regions. Then the first quartile (ADC L ) and last quartile (ADC H ) of ADC map were respectively overlaid on rCBV L maps. Finally, two intersections of ADC L -rCBV L and ADC H -rCBV L ROIs were obtained. Other regions within CE outside the two ADC-rCBV ROIs were taken as abnormal controls (CE control, CEC). Absolute volumes of ROIs were calculated in FSL [19]. Proportional volumes (%) of two ADC-rCBV ROIs were calculated as the ratio of the absolute volumes to CE volume.

MRSI voxel selection
Due to the difference in spatial resolution between MRSI and MRI, voxels from T2-weighted MRIs were projected onto MRSI using MATLAB, to evaluate the spectroscopic measure of the ADC-rCBV compartments identified on T2 space. The proportions of the voxels of the ADC-rCBV compartments occupying each MRSI voxel was calculated, and only the MRSI voxels that were completely within the delineated tumor were included in further analyses. Since the ADC-rCBV compartments potentially included multiple MRSI voxels, the proportions of T2-space voxels in the MRSI voxels were calculated and taken as the weight of each MRSI voxel. The sum weighted value was used as the final metabolic value, providing an objective method for MRSI voxel selection ( Supplementary Fig. 2).

DTI invasive phenotypes
We investigated DTI invasive phenotypes of 64 patients which overlap with a previously reported cohort and have been corre-lated to isocitrate dehydrogenase (IDH) mutation status [25]. Three invasive phenotypes ( Supplementary Fig. 3) were classified using previously described criteria [26] based on the decomposition of diffusion tensor into isotropic (p) and anisotropic components (q): (a) diffuse invasive phenotype; (b) localized invasive phenotype; and (c) minimal invasive phenotype.

Statistical analysis
All analyses were performed with RStudio v3.2.3. Continuous variables were tested with Welch Two Sample t-test. MRSI data or tumor volume, were compared with Wilcoxon's rank sum test or Kruskal-Wallis' rank sum test, as appropriate, using Benjamini-Hochberg's procedure for controlling the false discovery rate in multiple comparisons. Spearman's rank correlation was used to model the relation between the volume of two ADC-rCBV ROIs and the volume of CE and FLAIR ROIs. Survival was analyzed on patients who received standard treatment following Stupp's protocol. Kaplan-Meier's using log-rank test and Cox proportional hazards regression analyses were performed to evaluate patient survival. For Kaplan-Meier's analysis, the volumes of ROIs and MRSI variables were dichotomized using the 'surv_cutpoint' function in R Package ''survminer". Patients who were alive in last follow-up were censored. Multivariate Cox regression with forward and backward stepwise procedures was performed, accounting for relevant covariates, including IDH-1 mutation, MGMT promoter methylation status, sex, age, extent of resection and contrast-enhancing tumor volume. The forward procedure started from the model with one covariate. The backward procedure initiated from the model including all covariates. For each step, Akaike Information Criterion was used to evaluate the model performance. The final multivariate model was constructed using the covariates selected by the stepwise procedures. The hypothesis of no effect was rejected at a two-sided level of 0.05. Fig. 1. Illustration of the pipeline to identify two ADC-rCBV compartments. Both ADC and rCBV maps are co-registered to the T2 weighted images and tumor regions are segmented manually. Low perfusion regions are partitioned using a quartile threshold. Similarly, two ADC regions are partitioned using high and low ADC quartile thresholds respectively. The spatial overlap between the thresholded rCBV and ADC maps defined two compartments ADC H -rCBV L and ADC L -rCBV L . MR volumetric and metabolic analyses of both compartments are performed and interrogated in invasive phenotype and patient survival analysis models.

Multiparametric MRI identified two low perfusion compartments
The volumes of ROIs for patient subgroups are compared in Table 1. The interrater variability of the ROIs showed excellent agreement, with Dice scores of 0.85 ± 0.10 (CE) and 0.86 ± 0.10 (FLAIR) respectively. The ADC H -rCBV L compartment (volume 5.7 ± 4.6 cm 3 ) was generally larger than the ADC L -rCBV L compartment (volume 2.3 ± 2.2 cm 3 ) (P < 0.001). Completely resected tumors had smaller CE volume (P = 0.006) and smaller ADC H -rCBV L compartment (P = 0.002). Fig. 2(A, D) shows the two compartments for two cases.

Low perfusion compartments displayed abnormal metabolic signatures
The ADC H -rCBV L compartment showed a significantly higher lactate/creatine (Lac/Cr) ratio than NAWM (P < 0.001), and an increased ML9/Cr ratio compared to NAWM (P < 0.001). Similarly, the ADC L -rCBV L compartment displayed higher Lac/Cr ratio and ML9/Cr ratio than NAWM (both P < 0.001). Although not significant, the Lac/Cr and ML9/Cr ratios in the ADC H -rCBV L compartment were higher than the ADC L -rCBV L compartment (Supplementary Table 1). Fig. 2 shows the comparison of the metabolite levels of two compartments.

Low perfusion compartments exhibited diverse effects on tumor invasion
The contrast-enhancing (CE) tumor volume was significantly correlated with the Lac/Cr ratio in the ADC L -rCBV L (P = 0.018, rho = 0.34). Interestingly, the volume of tumor infiltration beyond contrast enhancement, which was delineated on FLAIR images and normalized by CE volume, showed a moderate positive correlation with the proportional volume of the ADC L -rCBV L compartment (P < 0.001, rho = 0.42) and a negative correlation with the proportional volume of the ADC H -rCBV L compartment (P < 0.001, rho = À0.32). The correlations of ROI volumes are demonstrated in Supplementary Fig. 4.

The ADC L -rCBV L compartment of minimally invasive tumors showed lower lactate
The minimally invasive phenotype displayed a lower proportional volume of ADC L -rCBV L compartment in CE ROI than the localized (P = 0.031) and diffuse phenotype (not significant), and a higher proportion of ADC H -rCBV L compartment than the localized (P = 0.024) and diffuse phenotype (not significant), suggesting the effects of the two low perfusion compartments to tumor invasiveness were different. Of note, the minimally invasive phenotype displayed lower Lac/Cr ratio compared to the localized (P = 0.027) and diffuse phenotype (P = 0.044), indicating that the ADC L -rCBV L compartment may have less acidic microenvironment in the minimally invasive tumors. A full comparison between the three invasive phenotypes can be found in Supplementary Table 2. Next, we included the volumes of two compartments and their Lac/Cr ratios into the survival models. The results using a stepwise procedure showed that higher volumes of the two compartments were associated with better PFS (ADC H -rCBV L : HR = 0.102,   Table 2 and the Kaplan-Meier curves using log-rank test are shown in Fig. 3.

Discussion
This study combined perfusion and diffusion parameters to quantify the low perfusion compartments that may be responsible for treatment resistance. The non-invasive approach using physiological imaging may potentially improve the commonly used weighted structural imaging.
The clinical values of the individual markers have been assessed previously. Among them, rCBV is reported to indicate IDH mutation status and associated with hypoxia-initiated angiogenesis [27]. Decreased ADC reflects restricted diffusivity, which is considered to represent higher tumor cellularity/cell packing [28] and associated with shorter survival [29]. Of note, although another meta-analysis also showed that the ADC value had an inverse correlation with cellularity in glioma, this correlation was not consistent in all tumor types [30]. Here we integrated multiparametric MRI to identify two low perfusion compartments. With similar low levels of perfusion, the restricted diffusivity in the ADC L -rCBV L compartment suggests this compartment may contain more microstructures, compared to the ADC H -rCBV L compartment. We measured the lactate, macromolecule and lipid levels at 0.9 ppm (ML9) in the spectra, as increased lactate indicates an acidic microenvironment, while ML9 is associated with proinflammatory microglial response [31]. The elevated ML9/Cr ratios may suggest both compartments displayed elevated inflammation response [31], potentially due to recruitment of inflammatory cells by necrotic tissue [32]. The positive correlation between tumor volume and lactate levels in the ADC L -rCBV L compartment could   indicate a higher lactate production as tumor grows. When evaluating the non-enhancing peritumoral regions, we found that tumors with larger infiltration area tended to have smaller ADC H -rCBV L and larger ADC L -rCBV L compartments, suggesting the latter might be more responsible for infiltration. This was supported by our findings that minimally invasive phenotypes displayed significantly lower lactate in the ADC L -rCBV L compartment. We further investigated the effects of two compartments on patient survivals. Interestingly, a higher Lac/Cr ratio in the two compartments was related to elevated hazard (HR > 1) while this ratio in other tumor regions showed a reduced hazard. This implies that the resistant phenotype may possibly reside in the two compartments. As the ADC L -rCBV L compartment was associated with larger tumor infiltration area, higher diffusion invasiveness, and the lactate level in this compartment had significant effect on both PFS and OS, this compartment may be associated with treatment resistance.
We found that higher volumes of both low perfusion compartments were associated with better survivals, while higher Lac/Cr ratios of these compartments were associated with worse survivals. These results suggested the higher proportion of the low perfusion compartments may indicate relatively lower tumor proliferation; the higher levels of lactate production in these compartments, however, may indicate a more aggressive phenotype. Consistent with our findings, previous studies showed that higher levels of lactate production were associated with radioresistance [33] and worse patient survival [34], potentially due to the antioxidative capacity of lactate.
Our findings have clinical significance. The identification of possibly resistant compartments could inform choice of treatment target. Our results show that higher acidic stress in ADC L -rCBV L compartment may lead to a more aggressive phenotype. Since adjuvant therapies may aggravate the microenvironmental stress, this finding suggests that more attention may be needed for patients with larger volumes of ADC L -rCBV L compartment.
There are limitations in our study. Firstly, due to the low spatial resolution of MRSI, the multivariate analysis was based on a subset of patients. The low spatial resolution of MRSI can also affect the comparison of the metabolic signatures of the two ADC-rCBV compartments. However, we have used a weighting method to include multiple MRSI voxels containing the two compartments, which may potentially help to reduce this bias caused by a single voxel containing both compartments. Secondly, the cut-off values defining the two compartments were based on the quartiles of rCBV and ADC distributions in individual lesion. We hypothesized that each individual tumor can be an independent ecological system in which the selective stress may arise from the disparities in sub-regional perfusion and diffusion. Compared to our method, a global absolute threshold across all patients may have the advantage of identifying consistent tumor habitats among patients. Several limitations, however, may still potentially exist when using the absolute threshold, even if we have normalized all the pixel values to the contralateral normal-appearing whiter matter. It may not address our hypothesis of intratumoral habitats with evolutionary disadvantages. More importantly, it could be significantly affected by the profound tumor heterogeneity and limited by the scanning setting used in this specific cohort. Thirdly, when the sequences of this study were designed, there was still no consensus regarding the use of pre-bolus of contrast agent. Therefore, a pre-bolus was not given in this study. However, to address this issue, we have used NordicICE to perform leakage correction across all patients using the same software setting. Lastly, although the imaging markers are validated histologically from other studies [35], a full biological validation can only be achieved with multi-region sampling of each tumor.
In conclusion, we showed that multiparametric habitat imaging could identify two low perfusion compartments, which may help improve the non-specific conventional imaging. The compartment demonstrating both low perfusion and restricted diffusion may indicate a habitat especially responsible for treatment resistance. This could provide crucial information for personalized treatment. As our analyses were based on clinically available imaging modalities, this approach could easily be implemented, and potentially extended to other system.