Introduction

Tumor heterogeneity poses a major clinical challenge in brain tumors and its enhanced assessment can be pivotal to understand treatment failure in high-grade gliomas (HGGs) (O’Connor et al., 2015; Stadlbauer et al., 2018b). Indeed, in HGGs, the regional diversity both at the level of DNA mutation and gene expression profiles may ultimately alter phenotype by influencing tumor cell metabolism, apoptosis, and angiogenesis processes (Louis et al., 2016; Komori et al., 2018). At the individual patient level, the heterogeneous biology of malignant gliomas challenges tools that rely on single biopsy criteria for making disease-wide assessments and poses a problem in planning effective subsequent treatments. Furthermore, tumor heterogeneity is reflected in differential and dynamical response during therapy, eventually causing the failure of previously effective alkylating drugs or targeted therapies, as resistant clones invariably emerge and proliferate (Sottoriva et al., 2013; Stadlbauer et al., 2018b).

Non-invasive imaging can provide multiple biomarkers with insight into the malignancy and biology of gliomas. Conventional magnetic resonance imaging (cMRI) is routinely employed in the diagnosis and clinical management of malignant gliomas: beside the tumor features described by cMRI sequences, advanced physiological MRI techniques such as diffusion MRI (dMRI) and perfusion-weighted imaging (PWI) add important structural, physiological and hemodynamic information to measure biological properties quantitatively and non-invasively and correlate with patients’ outcome (Castellano and Falini, 2016). Similarly, positron emission tomography (PET) imaging reflects fundamental metabolic patterns throughout the tumor. In particular, PET using tracers from nitroimidazoles family, such as 18F-labeled fluoroazomycinarabinoside (18F-FAZA) or 18F-labeled fluoromisonidazole (18F-FMISO), identifies areas of hypoxic tissue that have been shown as related to treatment resistance, thus negatively impacting on patient outcome and survival (Preibisch et al., 2017; Abdo et al., 2019).

Furthermore, with the application of advanced mathematical modeling, it is also possible to divide neoplasms into definite subregions comprehending clusters of voxels with comparable radiomics features. This method is defined as “habitat imaging” and is considered a way to characterize and measure the heterogeneity of the tumoral microenvironment (Gatenby et al., 2013). Habitat imaging has been currently mainly employed on cMRI sequences, although it could hypothetically be applied to any neuroimaging method, comprising advanced MRI techniques such as PWI and dMRI as well as PET imaging with different radiotracers (Zhou et al., 2014; Hu et al., 2015; Lee et al., 2015a,b; O’Connor et al., 2015; Cui et al., 2016; Sauwen et al., 2016; Dextraze et al., 2017; Kim and Gatenby, 2017; Zhou et al., 2017; Fathi Kazerooni et al., 2018; Fuster-Garcia et al., 2018; Ismail et al., 2018; Juan-Albarracin et al., 2018; Li et al., 2018a,b; John et al., 2019; Juan-Albarracin et al., 2019; Li et al., 2019; Stringfield et al., 2019; Wang et al., 2019; Wei et al., 2019; Zhang et al., 2019; Beig et al., 2020; Del Mar Alvarez-Torres et al., 2020; Verma et al., 2020; Wu et al., 2020; Kim et al., 2021; Park et al., 2021).

To date, the spatial relationship between areas of hypoxia, neoangiogenesis, and altered tissue microstructure and cellularity has not been completely elucidated in HGGs using advanced neuroimaging. The aim of the present study is to evaluate an innovative, PET and MRI approach for the assessment of hypoxia, perfusion, and tissue diffusion in HGG and derive a combined spatial habitat map for clustering of intra-tumor heterogeneity. To the best of our knowledge, this represents the first study combining the aforementioned parameters to define brain tumor habitats.

Materials and Methods Patients

In this study the PET and MRI datasets of twenty patients harboring brain lesions with MRI features suggestive for HGG were analyzed, previously recruited in a prospective study on hypoxia tracers at San Raffaele Scientific Institute from April 2016 to October 2017. The local research ethics committee approved the study, and informed consent was obtained from all patients. All the subjects underwent advanced MRI acquisition and 18F-FAZA PET/CT before surgery, except for two patients that were withdrawn from the study due to technical problems related to radiotracer synthesis. Another patient was eventually excluded from the analysis since the final histopathological diagnosis was brain metastases from lung cancer.

PET/CT

The procedure of 18F-FAZA radiotracer production followed the technique previously reported by Savi et al. (2017). Every recruited patient had an intravenous administration of about 372 ± 17 (minimum 340, maximum 407) MBq of 18F-FAZA. The PET/CT scan acquisition and analysis were conducted as previously reported by Mapelli and Picchio (2020) and Mapelli et al. (2021).

MRI Acquisition

Pre-operative MRI datasets were acquired (1–5 days prior to the surgical procedure) on a 1.5 T scanner (Philips Achieva–Philips Healthcare, Best, Netherlands) in 12/17 cases, and on a 3.0 T (Ingenia CX scanner, Philips Healthcare, Best, Netherlands) on 5/17 patients.

Imaging acquisition included the following sequences:

1. 3D-FLAIR (fluid-attenuated inversion recovery) and Contrast-enhanced (CE) 3D-T1 weighted images, for morphological evaluation and tumor segmentation (Figures 1A,B);

FIGURE 1

Figure 1. FLAIR-tumor and CE-tumor segmentations on different MRI sequences. (A) T1-weighted contrast-enhanced imaging: the CE-tumor volume is contoured in green, including the enhancing tissue and necrotic core only. (B) FLAIR images: FLAIR-tumor is contoured in orange, including both the CE-tumor and the peritumoral hyperintensity. (C) MD map elaborated from DTI acquisition. (D) rCBV map elaborated from DSC MR perfusion acquisition. (E) VP map elaborated from DCE MR perfusion acquisition (F) Ktrans map elaborated from DCE MR perfusion acquisition. (G) 18F-FAZA-PET uptake showing hypoxic areas. (H) Prototype of clusters resulting from the application of the spatial habitat imaging algorithm only to the contrast-enhanced tumor volume (CE-tumor); clusters obtained this way were considered too small and fragmented in order to be reliably correlated to histology, and therefore not suitable for the purpose of the present study. (I) Habitats obtained applying the spatial habitat imaging algorithm to the FLAIR-tumor volume (which included the contrast-enhanced tumor core along with the surrounding area of FLAIR hyperintensity) were able to divide the whole tumor in a more reproducible way, with larger clusters amenable to a reliable histological correlation with the bioptic samples obtained; for the aforementioned reasons we selected FLAIR-tumor as the final tumor mask to be investigated and all the results presented in our paper are based on its analysis.

2. Diffusion Tensor Imaging (DTI) sequence to detect microstructural features of the tumor tissue. DTI data were obtained using a single-shot echo planar sequence with diffusion gradients applied along 15 or 32 axes, using b-values of 0 and 1000 s/mm2. After correcting for distortions due to eddy currents and head motion, maps of Mean Diffusivity (MD) (Figure 1C) were derived from DTI acquisitions by fitting the diffusion tensor model within each voxel using the FMRIB Diffusion Toolbox (FDT tool, FMRIB Software Library [FSL] version 6.0.01);

3. Perfusion-weighted magnetic resonance imaging including dynamic contrast-enhanced (DCE) and dynamic susceptibility contrast-enhanced (DSC) MRI, to characterize tumor vascularization and neoangiogenic processes, according to a previously descripted protocol (Santarosa et al., 2016; Anzalone et al., 2018). PWI analysis was performed with Olea Sphere (v 3.0, Olea Medical Solutions, France) to obtain the parametric maps of volume transfer constant (Ktrans) and plasma volume (Vp), derived from DCE, and relative cerebral blood volume (rCBV), derived by DSC, as described in Conte et al. (2017) (Figures 1D–F). Pre-processing steps included automatic motion correction by a rigid-body registration, automatic spatial smoothing, and background segmentation.

The complete acquisition parameters of all the sequences are reported in the Supplementary Tables 1, 2.

Image Processing

For each patient, all the radiological studies and maps were imported, coregistered to the reference 3D-T1 CE sequence (two consecutive automatic rigid registration processes, with 6 degrees of freedom and 0.1 as percentage of samples, were performed after initial manual adjustment, then carefully checked for the eventual necessity of minor manual correction) and reformatted to the same matrix (voxel 1 × 1 × 1 mm) in the image processing software 3D Slicer2 (Fedorov et al., 2012), creating a single dataset for each patient. Two different tumor masks were segmented manually using 3D Slicer, being carefully checked repeatedly by the first author (MB) and an experienced neuroradiologist (AC): one including only the enhancing tissue and necrotic core on the post-gadolinium 3D-T1 images (CE-tumor, green contour in Figure 1A), the second including, additionally, the peritumoral infiltrative edema as identified on 3D-FLAIR images (FLAIR-tumor, orange contour in Figure 1B). Care was taken, while comparing data from multiple coregistered sequences (e.g., Susceptibility-weighted imaging, T2-weighted, PWI maps, MD maps, T1-weighted images with and without contrast enhancement, CT) not to include non-tumoral structures that could jeopardize subsequent analyses (e.g., subarachnoid cisterns, sulci, ventricular space, choroid plexus, ependyma, significant vascular structures, dural folds, calcifications).

The feasibility of deriving a combined spatial habitat map was then exploited by using the automatic Otsu thresholding algorithm method developed in Matlab2019a (The MathWorks, Inc., Natick, MA, United States). The Otsu algorithm was chosen for its simplistic approach in dividing ROI tumoral voxels into two clusters with high and low-intensity (Zhou et al., 2014, 2017; Napel et al., 2018; Stringfield et al., 2019), thus maintaining the clinical significance of each map subregion. The Otsu algorithm calculates a global threshold T from each of the grayscale tumoral ROI (perfusion, diffusion, and hypoxia, Figures 1C–G) by minimizing intra-class intensity variance, or equivalently, by maximizing inter-class variance as described in Otsu (1979).

For each subject the three tumoral ROIs were imported in Matlab2019a and voxels included inside the tumor ROI [excluding those outside the range of mean ± 3 standard deviations, as described in Collewet et al. (2004)] were clustered into two disjoint subregions with the automatic Otsu’s algorithm. Therefore, given three tumor-masked maps (D1, D2, and D3), they were clustered as follows:

D→1[D,1⁢LD]1⁢H,D→2[D,2⁢LD]2⁢H,andD→3[D,3⁢LD]3⁢H

D∩1⁢LD=1⁢H∅,D∩2⁢LD=2⁢H∅,andD∩3⁢LD=3⁢H∅

In this way, 2 clusters for each map (“high” and “low” expression of a parameter) were obtained. Then, spatial habitat maps were composed in each patient by defining as a different ‘habitat’ each possible combination of clusters (see Supplementary Figure 1). Thus, the cluster that groups together all high-intensity cluster sequences (HHH) is defined as:

Habitat(HHH):[D∩1⁢HD∩2⁢HD]3⁢H

The different habitats were eventually represented with different colors (Figures 1H,I) to ease visual identification, yield regions with different combinations of perfusion, diffusion, and hypoxia, supposedly reflecting different physiological microenvironments within the tumor.

Surgical Procedure and Histopathological Evaluation

All patients underwent the surgical procedure at the Neurosurgery Department of the IRCCS San Raffaele Hospital (Milan, Italy). Patients were considered for maximum safe resection or stereotactic needle biopsy according to tumor location/extension and patient characteristics (e.g., age, comorbidities, and performance status).

In the case of stereotactic biopsy, the procedure commenced with the fixation of the MRI-compatible Leksell stereotactic frame to the patient’s skull (Model G, Elekta, Stockholm, Sweden), under mild intravenous sedation (midazolam) and local anesthesia (lidocaine and mepivacaine). The patient then underwent the acquisition of an axial 3D-T1 MRI after gadolinium administration (Voxel 1 × 1 mm, slice thickness 1 mm, Matrix 256 × 256 mm) using the appropriate fiducial system on a Philips Achieva 1.5 Tesla MR scanner (Philips Medical Systems, Best, Netherlands). In most of the cases, a 3D T2-weighted series was also obtained. All the relevant imaging and parameter maps were then imported into the stereotactic planning software Leksell SurgiPlan® (Elekta Instruments AB, Stockholm, Sweden) and coregistered to the newly acquired stereotactic images.

Since the final tumor habitats were not available to the neurosurgeon pre-operatively, the designation of target sites for biopsy sites and the needle path were planned according to the best clinical practice in order to maximize patient safety while ensuring final diagnosis. During the procedure, a Sedan needle (2.5 mm diameter and 10 mm of bioptic window) was used to acquire a minimum of two bioptic samples. Considering the macroscopical aspect and size of the sample obtained and the procedure-related risks/complications, a median of 3 cylindrical tissue biopsies (median length 8 mm, range 4–11 mm) were obtained. The accuracy of the bioptic sampling was eventually verified in all cases through the co-registration of the pre-operative imaging and planned trajectory with the post-operative CT images inside the SurgiPlan software; the stereotactic coordinates of the exact, final sites of the biopsy were then imported in the 3D Slicer Dataset.

The tissue samples obtained were directly fixed in a 10% formalin solution and referred to the pathology department, where they were processed the same day or the following one if the procedures were performed in the late afternoon. The final diagnosis was established according to the 2016 World Health Organization (WHO) classification (Louis et al., 2016). A descriptive analysis of the microscopic morphology of each sample was performed by an expert neuropathologist (MC) attended by a neurosurgeon (MB), repeating the whole evaluation at three distant time points (blinded from previously reported values and patient’s case) in order to enhance the reproducibility of the results. Additionally, to address the issue of the length variability among different bioptic specimens and to maximize the correspondence with the target point set on the stereotactic planning system, only the central portion (extending 5 mm along the greatest axis of the specimen in its largest section available) was analyzed in each sample.

The following parameters were semi-quantitatively graded similarly to a previously reported work (Mikkelsen et al., 2018), despite the adoption of a higher number of classes dividing each feature in the present analysis: cellularity, the extent of necrosis, and the number of hyperplastic vessels counted. Cellularity was graded with a score from 1 (seemingly normal representation) to 5 (extremely high cellular density). Classes adopted for necrosis (expressed as a rate over the whole area considered for each specimen) were the following: 0 (0%), 1 (≤10%), 2 (11–25%), 3 (26–50%), 4 (51–75%), and 5 (76–100%).

In the 10/17 patients who underwent needle biopsy, the stereotactic coordinates of the sampling sites were imported in the 3D Slicer Dataset as a 5 mm diameter spherical volume, similarly to what was previously described by other authors (Jun et al., 2006; Barajas R. F. et al., 2012; Gates et al., 2019; Rotkopf et al., 2020), in order to take into account the maximum reported accuracy error in Euclidean distance reported in the literature (Roth et al., 2018). Each histopathological specimen, labeled according to the order of resection, was then analyzed in comparison to the corresponding voxels and habitats in the 3D Slicer Dataset.

Statistical Analysis

Statistical analysis was performed with IBM SPSS Statistics, version 23.0 (IBM Corp., Armonk, NY, United States), Prism 9 version 3.1 (GraphPad Software, LLC) and Matlab2019a (MathWorks, Natick, MA, United States). The threshold of statistical significance was defined as two-sided P = 0.05. Distribution normality was assessed with the Kolmogorov–Smirnov test.

The following variables were tested for possible associations: quantitative values across different habitats for each parametric map, CE-tumor volume, FLAIR-tumor volume, representation of the different habitats (absolute volume, rate over total), patient’s age, WHO tumor grading, MGMT methylation, 1p19q codeletion, IDH-1 status. The distribution of mean values of Vp, 18F-FAZA and MD parametric maps in each habitat was displayed by means of boxplots.

A Kruskal–Wallis test was conducted to determine if there were differences among different habitats in voxel-wise parametric values for each PET and MRI map. The percentage of habitats’ representation according to tumor grade and molecular status was assessed by boxplots’ visual inspection. Pairwise comparisons were performed using Dunn’s procedure with Bonferroni correction for multiple comparisons.

Results Patients and Lesion Characteristics

The study population (data summarized in Table 1) was composed of 5 female and 12 male patients, with ages ranging from 41 to 81 years old (mean 66 years).

TABLE 1

Table 1. Patient demographic, histopathologic, and survival data.

A single contrast-enhancing (CE) lesion was detected in 13/17 patients; 3/17 patients harbored multiple CE lesions, while in one patient no definite CE tumor core was confidently identifiable. The average cumulative CE-tumor volume (per patient) was 25.3 ± 24.4 cm3 (range 0–79.1 cm3). Mean cumulative FLAIR-tumor volume (per patient) was 94.2 ± 58.7 cm3 (range 10.2–185.4 cm3). The mean rate between CE-tumor and FLAIR-tumor volume was 0.30 ± 0.20 (range 0–0.7).

Maximum safe resection was performed in 7/17, while a stereotactic biopsy was obtained in 10/17 patients. Four out of individuals harbored an anaplastic astrocytoma (WHO grade III), one had an anaplastic oligodendroglioma (1p/19q co-deleted, WHO grade III) and 12/17 patients had grade IV glioblastoma, according to WHO 2016 classification (Louis et al., 2016). Three tumors retained an IDH-1 mutation. Histopathological and molecular data are summarized in Table 1.

Mean overall survival (Table 1) was 53.1 ± 67 weeks (median 21 weeks, range 8–222 weeks). Two patients were still alive at the last follow-up assessment (49 and 51 months, respectively). The remaining 15 patients died during follow-up due to causes directly (e.g., tumor progression increasing mass effect on cerebral structures) or indirectly (e.g., aspiration pneumonia) associated with the disease.

Map Selection

The ideal number of maps to be included simultaneously in the algorithm was thoroughly evaluated. The possible complementary information gathered by adding a higher number of maps had to be balanced with the exponential increase in the number of significantly smaller and fragmented clusters, resulting in an unreliable surgical targeting and histological correlation. For the aforementioned reasons, after examining the possible complementarity of the maps, “8” (from the combination of three different parametric maps) was identified as the ideal number of clusters for the purpose of the project. Among the available parametric perfusion maps, the intra-vascular plasma volume derived from DCE-MRI, that quantifies tumor neoangiogenesis similarly to rCBV (Alcaide-Leon et al., 2015), was selected as the most appropriate for the purpose of the study due to the lower sensitivity to susceptibility artifacts and to the presence of intralesional hemorrhagic foci, as well as for the better capability to discriminate surrounding macro-vessel from true intralesional neoangiogenic hotspots on Vp parametric maps, especially in cortical tumors growing in proximity to large vessels (Santarosa et al., 2016). Furthermore, DCE-derived Vp maps have higher spatial resolution compared to DSC-derived rCBV maps and yield an absolute quantification of perfusion parameters (Anzalone et al., 2018). Vp also allows a more informative representation of the FLAIR-tumor volume outside of the CE region compared to Ktrans (Nguyen et al., 2016; Zhang et al., 2017; Anzalone et al., 2018).

Positron emission tomography (PET) with the 18F-flouroazomycin arabinoside (18F-FAZA) radiotracer has been previously used for the non-invasive assessment of tumor hypoxia (Mapelli et al., 2021), and was chosen due to its high perfusion, high biological half-life time and fast clearance from blood, with subsequent improved tumor-to-background ratio (Mapelli and Picchio, 2020).

Among the DTI-derived maps, MD was selected as one of the simplest and most validated parameters in radiomics to distinguish tumor grade and cellular proliferation in brain gliomas (Miloushev et al., 2015). MD quantifies water molecular mobility, averaged over all the obtained directions, within tissues impeded by cell membranes and tortuosity of the extracellular space (Maier et al., 2010). This parameter is equivalent to the apparent diffusion coefficient (ADC) generated from DWI acquisitions and shows a negative correlation with tissue cellularity, as previously reported in literature (Ellingson et al., 2010; Maier et al., 2010; Chen et al., 2013; Sanvito et al., 2021).

Habitats’ Analysis

Eight habitats were generated from cluster intersection by combining all possible HIGH (H) and LOW (L) intensity regions of each map (Vp + 18F-FAZA + MD) and the final tumor habitat map was composed (Figure 1).

In the preliminary analysis of our experiment, we investigated the application of the same clustering algorithm to either the CE- or FLAIR-tumor, since both of them have been considered for habitat imaging studies previously published, and we assessed the obtained habitats. Clusters obtained from CE-tumor volume alone were considered too small and fragmented in order to be reliably correlated to histology, and therefore they were not deemed suitable for the purpose of the present study (Figure 1H). On the contrary, habitats obtained from FLAIR-tumor volume (which included the CE region along with the surrounding FLAIR hyperintensity) were able to divide the whole tumor in a more reproducible way, with a larger cluster amenable to a reliable histological correlation with the bioptic samples obtained (Figure 1I). For these reasons, we selected FLAIR-tumor as the final tumor mask to be investigated and all the results presented in our paper are based on its analysis. The radiological characteristics of the 8 clusters, and their volumetric representation, are summarized in Table 2 and depicted in Figure 2. Median voxel values among different maps (Vp, 18F-FAZA, and MD) and Habitats (1–8) are reported in Table 3. Boxplots of the quantitative distribution of mean values of Vp, 18F-FAZA and MD parametric maps in each environment are shown in Figure 3.

TABLE 2

Table 2. Habitats’ characteristics and volumetric representation.

FIGURE 2

Figure 2. Habitats volumetric representation. The relative volumetric representation of each habitat (expressed as a percentage of the whole FLAIR-tumor volume) is reported for each case included in the series.

TABLE 3

Table 3. Plasma volume (Vp), 18F-labeled fluoroazomycinarabinoside (18F-FAZA), and mean diffusivity (MD) voxel median values among different habitats.

FIGURE 3

Figure 3. Perfusion (Vp), 18F-FAZA, and Diffusivity (MD) quantitative values within each habitat for the whole group of patients. Boxplots shows the distribution of mean values of Vp, 18F-FAZA, and MD parametric maps in each environment.

Figure 4 depicts habitats distribution in four different cases of the present series. Figures 5, 6 present two illustrative cases showing the consistency between the functional spatial tumor habitats and the microscopic morphology of the corresponding tissue. The microscopic semi-quantitative analysis (cellularity, extent of necrosis, number of hyperplastic vessels) of each sample and the corresponding Habitats on imaging are reported in Table 4. Table 5 and Figures 7, 8 report habitats’ prevalence according to grade and molecular features. Habitats’ descriptions and characteristics are analyzed below.

FIGURE 4

Figure 4. Four different illustrative cases (A–D) showing habitats distribution and reproducibility of spatial habitat imaging maps. A table listing the key colors (1–red, 2–orange, 3–gold, 4–yellow, 5–aquamarine, 6–light blue, 7–dark blue, and 8–violet) assigned to each habitat is reported in Figure 5. The yellow arrows shows point areas of microhemorragic foci resulting in a MD signal restriction.

FIGURE 5

Figure 5. Illustrative case of a 57-year-old female patient (case 16 in Table 1) who underwent stereotactic biopsy. (A) The trajectory path (yellow line) and the four consecutive sampling sites (enumerated in the yellow circles) are displayed in axial formatted pre-operative T1-CE, FLAIR, MD, 18F-FAZA, Vp, and post-operative CT images; the spatial distribution of the identified habitats and the assigned key colors are also displayed (1–red, 2–orange, 3–gold, 4–yellow, 5–aquamarine, 6–light blue, 7–dark blue, and 8–violet). (B) Microscopic slide [20× magnification, Hematoxylin and Eosin (H&E)] from the first sampling (corresponding to Habitat 7–Low Perfusion/Low Hypoxia/High Diffusivity and Habitat 8–Low Perfusion/Low Hypoxia/Low Diffusivity) showing a cortical/subcortical area with normal tissue architecture, presence of normal vessel and no significant increase in cellularity, but rather mild reactive gliosis. (C) Microscopic slide (20× magnification, H&E) from the second sampling (Habitat 1–High Perfusion/High Hypoxia/High Diffusivity) showing a clearly pathological area of white matter with increased cellularity, disruption of the physiological tissue architecture, and abnormal vessels. (D) Microscopic slide (20× magnification, H&E) from the third sampling (Habitat 2–High Perfusion/High Hypoxia/Low Diffusivity) showing neoplastic tissue with cytoarchitecture disruption, elevated cellularity with the presence of gemistocytic astrocytes, nuclear atypia, abundant eosinophilic cytoplasm, abnormal vascularization. (E) Microscopic slide (20× magnification, H&E) from the fourth sampling (Habitat 5–Low Perfusion/High Hypoxia/High Diffusivity and Habitat 7–Low Perfusion/Low Hypoxia/High Diffusivity) showing tissue with extensive necrosis.

FIGURE 6

Figure 6. Illustrative case of a 58-year-old male patient (Case 2 in Table 1) who underwent stereotactic biopsy. (A) The trajectory path (yellow line) and the three consecutive sampling sites (enumerated in the yellow circles) are displayed in axial formatted pre-operative T1-CE, FLAIR, MD, 18F-FAZA, Vp, and post-operative CT images; habitats were identified using the same key-color described in Figure 5. (B) Microscopic slide [20× magnification, Hematoxylin and Eosin (H&E)] from the first sampling (Habitat 7–Low Perfusion/Low Hypoxia/High Diffusivity and Habitat 8–Low Perfusion/Low Hypoxia/Low Diffusivity) showing a seemingly normal white matter architecture. (C) Microscopic slide (20× magnification, H&E) from the second sampling (Habitat 1–High Perfusion/High Hypoxia/High Diffusivity and Habitat 2–High Perfusion/High Hypoxia/Low Diffusivity) showing evident hypercellularity with nuclear atypias and marked vascular proliferation. (D) Microscopic slide (20× magnification, H&E) from the third sampling (Habitat 5–Low Perfusion/High Hypoxia/High Diffusivity and Habitat 6–Low Perfusion/High Hypoxia/Low Diffusivity) areas of vital tumoral tissue (on the upper part) mixed with areas of necrosis (on the lower part).

TABLE 4

Table 4. Semi-quantitative histopathologic analysis for each bioptic sample obtained.

TABLE 5

Table 5. Mean and median representation of the eight different habitats according to IDH-1 status and WHO 2016 grade.

FIGURE 7

Figure 7. Habitats’ representation (% of total tumor volume) according to WHO 2016 grade. For each habitat, the boxplot on the left (with filled color) represents WHO grade III tumors, while the boxplot on the right (not filled with color) represents WHO grade IV tumors. WHO grade IV gliomas had a higher (although not statistically significant) representation of Habitat 1 and 2 (considered the most aggressive clusters from a metabolic point of view) and Habitat 5 (mostly associated with necrotic areas). The most aggressive clusters are surrounded by a gray box with a dashed line.

FIGURE 8

Figure 8. Habitats’ representation (% of total tumor volume) according to IDH-1 status. For each habitat, the boxplot on the left (with filled color) represents IDH-1 Wild-type tumors, while the boxplot on the right (not filled with color) represents IDH-1 mutated tumors. Main differences were found in the representation of Habitat 7 (possibly embodying pure vasogenic edema) and 8 (possibly representing infiltrative edema). The former was significantly (P = 0.05) higher in IDH wild-type tumors, while the latter was higher (not significantly) in IDH-mutated lesions; these latter two habitats were surrounded by a gray box with a dashed line to highlight the differences.

Habitat 1, characterized by high perfusion (the highest Vp values among all clusters, Figure 3), high hypoxia, and high diffusivity, corresponded to areas of mild hyperintensity on T2/FLAIR sequences and intense enhancement on T1–CE sequences; Habitat 1 voxels were localized inside of CE-tumor segmentation in 92% of the cases overall. No areas of necrosis were identifiable on cMRI sequences inside this cluster. This cluster was often located at the periphery of the CE-tumor, leaving other aggressive clusters more internally in the tumor core. Habitat 2 was denoted by elevated perfusion, hypoxia, and restricted diffusion, representing, alone, about half of CE-tumor volume overall. On cMRI it was denoted by a mild hyperintensity on T2/FLAIR sequences and no visible necrosis. Similarly to Habitat 1, this cluster was localized inside of CE-tumor segmentation in 89% of the cases overall, usually in close proximity to the previous one; however, when compared to the first cluster, Habitat 2 was much more represented (12.1% vs. 4.5%), harboring a less intense contrast enhancement on T1-CE sequences, lower perfusion values, more reduced diffusivity (possibly reflecting higher cellularity) and, generally (in 12/17 patients), higher hypoxia (18F-FAZA-uptake). At histopathology these first two habitats were largely associated with a higher number of hyperplastic vessels and cellularity, whereas a very low rate of necrosis. Moreover, their relative prevalence was higher among WHO grade IV gliomas.

Habitats 3 and 4, both sharing low hypoxia and high perfusion, while a different degree of diffusivity, were extremely fragmented and poorly represented among all cases (4.5 and 5.2% of the tumor volume, respectively). They were mostly grouped together, often located either at the interface between tumor CE core and the surrounding area of FLAIR-hyperintensity with no CE, adjacent to large vessels or close to the cerebral cortex. Only about 25% of the voxels of these clusters were located inside of CE-tumor segmentation.

Habitat 5 was characterized by high hypoxia, low perfusion, and increased diffusivity (possibly denoting low cellularity). It was located inside of CE-tumor segmentation in 93% of cases, characterized by necrotic areas at both cMRI and histopathology. This was the least prevalent habitat, with an average representation of 4.2% (range 0.3–12.4%) among the considered tumors; however, it appeared more represented in WHO grade IV (4%) compared to WHO grade III tumors (2.1%).

Habitat 6 (low perfusion, reduced diffusivity, and high hypoxia) had a relatively high representation compared to the aforementioned clusters (average 7.1%, range 0.7–23%), being mostly localized inside of CE-tumor segmentation (75% of the cases). It was predominantly identified as areas of less intense CE inside the tumor core, being closely related to Habitat 5 (located more internally) and Habitat 2 (located more externally), as well perceivable in Figures 4C,D. Conversely, in different cases, this cluster was identified as areas of faint contrast-enhancement outside the main tumor cores, as one could appreciate in Figures 4A,B. At histopathology this habitat showed features of intermediate aggressiveness.

The remaining two clusters had the largest volumetric representation (each one representing over 30% of the total FLAIR-tumor volume) and were both found in areas with no CE in roughly 95% of the cases. Habitat 7 was, possibly, the least aggressive one inside the tumor mask due to its low values in all the considered parameters; it was characterized by the highest absolute MD values among all clusters (possibly reflecting the lowest cellularity) and it was located in areas of markedly hyperintense FLAIR signal (Figures 4B,D). The distribution of Habitat 7 was significantly higher (mean 33.3% vs. 18.8%, median 34.7% vs. 15.4% P = 0.05, independent sampled Mann-Whitney U test) in IDH-1 wild-type tumor compared to IDH-1 mutated (see Table 5 and Figure 8). Habitat 8 was characterized, on the contrary, by reduced diffusivity; moreover, perfusion values and hypoxia values tended to be higher than those of Habitat 7 (Table 3). Habitat 8 was also identified in the setting of hemorrhagic foci or necrotic areas with micro-hemorrhagic components (as identified on susceptibility weighted imaging, SWI) or viscous mucinous components, due to low MD values resulting from extravascular blood components or liquefactive necrosis (Kang et al., 2001; Chen et al., 2018; Ono et al., 2018) (yellow arrow on Figures 4A,B). Additionally, as easily noticeable in Figures 4C,D, the localization of Habitat 8 in cortical areas distant from the tumor core is likely related to the lower MD values of the cerebral cortex, compared to underlying edematous white matter. When considering bioptic samples obtained outside the CE-tumor these last two habitats were characterized by histological features of low cellularity and no signs of necrosis or neovascularization. Conversely, they showed extensive necrosis with no viable tissue, in those cases biopsied inside the CE-tumor.

Discussion

In this study, quantitative radiomic parameters from PET and MR images were derived for the assessment of hypoxia, perfusion, and tissue diffusion in a cohort of high-grade glioma patients, and a pipeline for a voxel-wise analysis to combine them in a number of discrete spatial tumor habitats was implemented for the clustering of highly hypoxic areas as well as areas of similar perfusion and diffusion characteristics. By this approach, eight regions were consistently highlighted within the tumor volume, reflecting the intra-tumor heterogeneity of physiological and histopathological microenvironments. Spatial tumor habitats demonstrated high reproducibility across patients and correlation with disease-specific histopathological and cMRI features. To the best of our knowledge, this represents the first study combining the aforementioned parameters to define tumor habitats.

Combined PET and MRI Hypoxia, Perfusion, and Diffusion Spatial Habitats’ Features

Of the eight habitats identified, we might speculate that the first two could identify the “vital” core of the tumor, with the second cluster being, hypothetically, the most aggressive component. Interestingly, Habitat 1, characterized by the highest Vp values (Figure 3), was often located more peripherally than Habitat 2, which is in line with the fact that perfusion values tend to be higher at the periphery of the tumor core (Ponte et al., 2017). These two habitats also showed the highest hypoxia values, data in accordance with previously reported studies (Ponte et al., 2017). The fact that the tumor “vital core” represents, at the same time, the most perfused and most hypoxic area, although it might seem paradoxical, is likely related to the functionally and structurally abnormal vascularity that characterizes HGGs, unable to provide sufficient perfusion and oxygen supply (Ponte et al., 2017; Rosinska and Gavard, 2021). The higher hyperintensity on T1-weighted CE sequences noticed inside Habitat 1 could be related to higher perfusion and/or vascular permeability, as confirmed by the significantly higher Vp, but also rCBV and Ktrans values noted in this cluster. The higher MD values of Habitat 1, when compared to Habitat 2, could imply lower tissue cellularity but also the presence of stronger vasogenic edema due to blood-brain barrier leakage. As somehow expected, the distribution of these habitats, together with the other most “aggressive” biologic cluster representing areas of tumor necrosis (Habitat 5), was higher in WHO grade IV tumors compared to WHO III.

The poor representation of Habitat 3 and 4 might be related to the low prevalence of well-oxygenated areas inside the tumor core. Since a certain degree of hypoxia is required to trigger neoangiogenesis (and thus increased perfusion), it is unlikely to find highly perfused but scarcely hypoxic clusters within HGGs (Huang et al., 2016; Ponte et al., 2017).

Habitat 5 represents necrotic areas still including some degree of vital tissue, since the 18F-FAZA-PET tracer needs hypoxic, yet vital, cells in order to be reduced to reactive oxygen radicals and remain locally trapped (Halmos et al., 2014; Mapelli and Picchio, 2020).

Habitat 6, when colocalizing with areas of faint contrast-enhancement outside the tumor “vital” and necrotic cores, was interpreted by the authors as new tumoral seeds progressing toward higher malignancy, since it harbored reduced diffusivity and increased hypoxia, but not yet a significantly augmented vascularization. On the contrary, when Habitat 6 was identified inside of the tumor core, it was interpreted as regions where tumor growth exceeded neovascularization and that were close to evolving into necrosis.

Finally, we suggested that Habitat 8, due to its reduced diffusivity, might identify areas of tumor infiltration beyond the CE vital core, while Habitat 7 areas of vasogenic edema (with no, or less prominent, tumor infiltration). Interestingly, in our series, the relative volume of pure edema (Habitat 7) was significantly higher in IDH-1 wild-type tumors, compared to IDH-1 mutated lesions, similarly to other recent studies (Dubinski et al., 2021).

Considering the known models of tumor progression (Holash et al., 1999; Brat et al., 2002; Hardee and Zagzag, 2012; Stadlbauer et al., 2018b), we may hypothesize a possible evolution over the time of the identified PET/MRI habitats: areas of tumor infiltration in the white matter (Habitat 8) might eventually develop some degree of hypoxia (Habitat 6), which, in turn, trigger neovascularization possibly through tumor cells secretion of mediators such as HIF-1α and VEGF (Stadlbauer et al., 2018b) and increased perfusion values, becoming a tumor “vital core”: Habitat 1 and (perhaps subsequently) Habitat 2. In due course, the tumor growth in the vital core will overwhelm its abnormal vascularization, leading toward areas in a pre-necrotic state (Habitat 6); as the necrosis progresses, cellularity will decrease (Habitat 5) eventually turning into liquefactive necrosis, with no remaining vital cells and showing, consequently, a concomitant drop in hypoxia signaling and 18F-FAZA-tracer uptake (Habitat 7). The aforementioned “habitat progression” is easily perceivable in Figures 4C,D, where a somehow concentric distribution of the clusters is identifiable.

Potential Clinical Applications

Advanced neuroimaging techniques such as PWI, dMRI, or PET with different radiotracers are usually evaluated independently in the clinics or are rather combined into an average parameter for the entire region of interest, discarding important spatial information (O’Connor et al., 2015; Napel et al., 2018). The simultaneous combination of diverse imaging modalities into a single habitat map leverages the information of the individual imaging datasets, possibly enhancing the ability to capture subtle physiological differences within tumor tissue and better revealing the complexity of spatial and temporal heterogeneity, both at initial diagnosis and during follow-up (Yip et al., 2017; Napel et al., 2018).

Furthermore, the heterogeneous biology of HGGs challenges the reliability of limited needle biopsy sampling for making disease-wide characterization. Nevertheless, needle biopsies are utilized, to date, in about a third of cases to define HGG diagnosis. The rate of non-diagnostic (to be repeated) stereotactic biopsies is 7–10% (Lu et al., 2015; Sciortino et al., 2019), and complete diagnostic accuracy has been reported as low as 33–64% of cases (compared to pathologic analysis obtained with the subsequent gross-total resection of the same lesion) (Jackson et al., 2001; Muragaki et al., 2008; Reithmeier et al., 2013). Integrating spatial tumor habitats into a bioptic planning system might better guide the neurosurgeon in sampling the most relevant functional regions of the tumor, especially by virtue of the much easier perception/understanding of the colored maps of a limited number of habitats, compared to the thorough expertise in advanced neuroimaging needed to interpret multiple coregistered parametric maps elaborated from PWI, dMRI, or PET radiotracers (as depicted in the illustrative cases of Figures 5, 6). This might be particularly valuable in those gliomas showing poor or no enhancement on T1-CE sequences, since about one-third of non-enhancing gliomas contain aggressive areas that, if correctly included in the pathologic analysis, may lead to classify the lesion into a higher-grade tumor, with consequently remarkable differences on treatment management and prognosis. In these cases, the few samples achievable via needle biopsy might easily underestimate the actual tumor grade if not adequately guided to the most aggressive areas.

We believe that further exploration of hypoxia-, perfusion-, and diffusion-derived habitats might provide precious, holistic information about intra-tumor heterogeneity, unveiling regional diversity both at the level of DNA mutation and gene expression profiles, helping in prognosis stratification and, perhaps, in future, treatment selection. The precise identification of the areas with the highest perfusion, microstructural alteration and cellularity, and hypoxia inside and outside the tumor contrast enhancement component could represent an enormous aid in surgical and radiotherapy planning. In the former case, the identification of an aggressive habitat might induce the surgeon to seek its maximum safe resection rate even if not showing evident contrast enhancement. In the latter case, radiation dose painting with boosts on the most aggressive and resistant clusters represents one of the most appealing perspectives for the next future (Castellano et al., 2021). Finally, the identified habitats could represent a precious tool in monitoring the dynamical spatio-temporal modification in glioma heterogeneity during therapies and follow-up and could identify more clearly areas of tumor response or progression.

Limitations and Future Work

One main drawback of the present work is the limited number of cases analyzed and their heterogeneity (e.g., tumor histology, grading, molecular characteristics, volume, location, patient’s age, comorbidities, type of surgery) which have likely hampered the achievement of conclusive associations between the variables involved (parametric map values, habitats, histopathological data, and patient outcome/survival).

In the subgroup of 10 patients undergoing a stereotactic biopsy, we were able to identify a qualitative correlation between the expected microenvironment of the different habitats and the actual histopathologic morphology, although only a preliminary analysis was possible in this small and heterogenous series and no definitive, statistically significant, conclusion could be drawn. The designation of biopsy target spots and path was planned in order to maximize patient safety while ensuring final diagnosis, according to the best clinical practice, by means of the minimum amount of sampling deemed necessary, due to ethical concerns. This, along with the fact that the spatial tumor habitats were not available to the neurosurgeon pre-operatively, prevented the acquisition of a sufficient variety of sampled habitats in the same lesion and frequently ended in sampling over multiple clusters at once, therefore limiting the possibility of direct associations. A thorough evaluation of prospectively and systematically collected data is clearly necessary to correlate the identified clusters with histopathological, immunohistochemical, and molecular biomarkers.

Multiparametric analysis on a per-voxel basis to generate habitats maps remains challenging, because this requires spatial registration of all imaging sequences, which can be hampered by the resolution differences and deformations that can occur between different acquisition modalities (Winfield et al., 2016; Napel et al., 2018). In the future, the use of a fully hybrid PET/MR scanner could further enhance the quality and spatial discrimination of the data obtained from dMRI and PWI while optimizing scanning time and image coregistration. Moreover, considerable processing and analyses were required initially for the generation of the final spatial tumor habitats, which demanded intensive labor and both clinical imaging scientist and informatics/engineering expertise. This requirement could be particularly critical in clinical practice when a short time in-between imaging and surgery is available. However, we believe that this could be overcome in the future via dedicated flow-charts with automated algorithms, providing the benefit of a much easier data interpretation for the neurosurgeon in charge of planning the procedure.

A limitation of the applied algorithm was the possible colocalization of the same radiological habitat in different tumor microenvironments. This mainly derived from the fact that MD, despite it has been proved to be inversely correlated to tissue cellularity, also depends on other parameters such as free water, perfusion/ischemia, temperature, and the presence of blood or mucinous components (Maier et al., 2010; Chen et al., 2018; Ono et al., 2018). For instance, Habitat 8 (which was thought to mainly represent areas of increased cellularity due to tumor infiltration of the surrounding white matter) was also found in proximity to microhemorrhagic foci, liquefactive necrosis including viscous mucinous components or cerebral vessels and cortex. On the other hand, perfusion values may be partly influenced by the presence of macrovessels passing through the analyzed tissue. A possible solution to the issues above could be leaving a larger margin from the aforementioned regions while segmenting tumor mask, adding/modifying the parametric maps utilized in the procedure, or identifying different cut-offs or algorithms for the clustering analysis.

The adopted Otsu algorithm is an easy-to-apply method able to find the optimal binary threshold for a single parametric map and outputting two disjoint regions, with high and low-intensity voxels. This allows combining different map subregions to obtain clusters of possible clinical significance. However, this method is limited in capturing all of the variety of tissues present in the tumoral ROI due to the net cut-off of the binary threshold, and it inherently assumes that all maps contain only two tissue classes, neglecting the fact that some maps could be sensitive to multiple tissue microenvironments. Moreover, this technique quickly leads to challenges when trying to combine multiple maps: the addition of more than 3 maps to the computation was found inadequate of the purpose of the study, due to the exponential increase of clusters number, making them excessively small and fragmented within the tumor ROI to offer a reliable correlation with bioptic specimens. The adoption of more refined data-driven clustering methods such as k-means, which allow to consider multiple parametric maps at the same time while preventing over-segmentation would have, perhaps, better addressed this particular issue. However, for the purposes of the present work, the generated clusters needed to be consistent in number and quantitative features among the different cases, in order to allow a possible biological interpretation of the individual habitats.

The alternative adoption of novel advanced MR techniques capable of assessing tissue hypoxia could represent another extremely appealing eventuality, since it would avoid the need for PET hypoxia radiotracers (such as 18F-FMISO or 18F-FAZA) that can be difficult to be gathered from specialized manufacturers or produced locally, expensive, time-consuming, and still problematic to be implemented in common practice due to the lack of a thorough clinical validation and institutional authorization. Indeed, in more recent years, resting-state blood oxygen level-dependent (BOLD) acquisitions have been successfully employed in the detection of areas of neurovascular uncoupling and hypoxia inside of tumoral tissue (Stadlbauer et al., 2011, 2018a,b; Pillai and Zaca, 2012; Englander et al., 2018). These maps could be easily acquired in the same MRI protocol together with morphological sequences, dMRI, and PWI, and can be eventually integrated into our algorithm to generate spatial tumor imaging habitats.

This study represents a preliminary analysis, as part of a feasibility study: additional future work will involve performing a systematic correlation of spatial habitat analysis with patients’ outcome, survival and recurrence, in a larger, prospective clinical series. Moreover, it would be interesting to follow-up patients over the course of the disease with the tumor habitat maps in order to monitor spatio-temporal changes and be able to validate the speculated spatio-temporal evolution of the PET/MRI tumor habitats.

Conclusion

This study shows the feasibility of an innovative PET and MRI approach for the assessment of hypoxia, perfusion, and diffusion in HGGs to derive combined habitats for the clustering of intra-tumor heterogeneity. The high consistence and reproducibility of the proposed spatial habitat imaging approach, as well as the correlation of the identified habitats with disease-specific histopathological features suggest that this approach can be valuable to decode metabolic, structural, and physiological patterns at same moment and non-invasively, in order to identify different areas within heterogeneous tumors. A systematic, prospective validation of this combined PET and MRI hypoxia, perfusion, and diffusion spatial tumor habitat map in a larger cohort of malignant glioma patients is warranted to define its impact on individualized therapy planning, prediction of prognosis and follow-up.

Data Availability Statement

All data generated and analysed during the current study are available from the corresponding author on reasonable request.

Ethics Statement

The studies involving human participants were reviewed and approved by the Institutional Ethics Committee of Ospedale San Raffaele, Milan, Italy. Patients provided their written informed cons