Glioblastoma (GBM) is a diffuse and highly invasive brain tumor of astrocytic origin. Despite recent advances in medical imaging diagnostics, radiation treatment, and chemotherapy, gliomas remain the most lethal cancer of the central nervous system (1). The median survival time is 14.6 months with a 5-year survival rate of 5.1% (2, 3).
GBMs are pathologically defined by pseudopalisading necrosis (4) and endothelial proliferation (5) and are highly heterogeneous with areas of differing cellularity, vascularity, and necrosis. This heterogeneity is further influenced by differing cellular microenvironments (6–8). Recent studies have hypothesized that these microenvironments create unique radiographic signatures on magnetic resonance (MR) imaging that may indicate information of clinical importance (8–10).
Magnetic resonance imaging (MRI) is commonly used to diagnose GBM, monitor tumor progression, and assess response to therapy (11, 12). Multiparametric MRI (MP-MRI) combines different MR contrast sequences to provide additional complementary structural information, potentially illuminating hidden characteristics and offering insight into both tumor and normal tissues. Recently, techniques that use advanced image processing to extract and analyze quantitative imaging features have been used in combination with clinical variables. These techniques are broadly defined as radiomics (9, 10, 13–15). Features quantified include metrics such as intensity distribution, spatial relationships, textural heterogeneity, and shape descriptors, among many other characteristics (16–18). The resulting matrices of feature data can then be mined along with a clinical variable of interest, such as a genetic profile, to determine previously hidden patterns (9, 19, 20).
The goal of this study is to develop an intuitive methodology for extracting voxel-wise radiomic profiles (RPs) from clinical imaging to describe distinct intensity characteristics from a given combination of MRI contrasts. We hypothesized that heterogeneous RPs represent underlying tumor characteristics, and thereby contain prognostic significance. This retrospective study looks at imaging from patients with GBM before treatment and derives a model of risk stratification for patient prognosis based on volumetrically thresholded RPs.
Patient Population and Image Acquisition
In total, 81 patients with pathologically confirmed primary GBM with T1-weighted pre- and postcontrast images, diffusion-weighted imaging, and T2-weighted fluid-attenuated inversion recovery (FLAIR) imaging before therapy were included in this study. Overall survival (OS) was corrected to account for contrast-enhancing tumor size, which is known to be associated with OS. Imaging was gathered before surgery using one of our institutional MRI scanners, including 1.5 T and 3 T GE (General Electric Health, Waukesha, Wisconsin) and Siemens magnets (Siemens Healthcare, Erlangen, Germany). The following are example scan parameters (all given in the format repetition time/echo time) at 1.5 T: T1 spin-echo sequence, 666/14 milliseconds; contrast-enhanced T1 acquired with gadolinium, 666/14 milliseconds; apparent diffusion coefficient, calculated from diffusion-weighted images acquired with an echo planar/spin echo sequence, 10 000/90.7 milliseconds; and FLAIR, acquired with an inversion recovery sequence, 10 002/151.8 milliseconds and TI of 2200 milliseconds. All images were acquired with submillimeter in-plane resolution.
Images were coregistered to the T1 image using FSL's FLIRT command (FMRIB Software Library, Oxford) (21–23). The FLAIR hyperintensity and contrast-enhancing lesion were semiautomatically segmented using a thresholding technique followed by manual correction of misclassified voxels (24). To correct for subtle intensity variance, images were smoothed with a 2-mm full width at half maximum Gaussian filter. Each image was then intensity-normalized by dividing voxel intensity by the standard deviation of the whole brain (25, 26). A total region of interest (ROI) mask was created using the union of the contrast-enhancing lesion and the FLAIR hyperintensity ROIs. Each contrast was then masked to only include abnormal regions. FSL's FAST command (FMRIB Software Library, Oxford) (27) was used to create a 3 tissue-class segmentation of each of the 4 images, ranking voxel intensities by low, medium, and high. Each was then coded with a value of 1, 2, and 3, respectively.
Images were coded based on the RP on a voxel-wise basis; for ease of interpretation, a 4-digit code was assigned to each voxel representing the intensity-based segmentation from each of the 4 images. The digit order chosen was T1, ADC, T1+C, and FLAIR. Codes ranged from 1111, representing dark voxels on all images, to 3333, representing all bright voxels. The left panel of Figure 1 shows the 4 clinical images used to generate RPs. The ROI used for all 4 images is the T1 enhancement combined with the FLAIR hyperintensity. These images are segmented individually within the ROI and given values of 1–3 based on intensity and neighboring pixel information. These 4 images are combined into one, with each voxel containing the segmentation value from all 4 images.
An RP is defined as all voxels within an image that contain the same 4-digit code. The resulting map contains 81 (34) potential RPs. Profiles were evaluated within the following 4 ROIs:
We performed an exploratory analysis to first determine which RPs were correlated with OS. The optimal volumetric threshold was then calculated for each significant profile using a log-rank Kaplan–Meier survival analysis. High- and low-volume groups were required to have at least 10 patients. Because of the high number of statistical tests performed, a strict P value of <.0005 was considered significant.
To generate a combined indicator score, each patient was given a score from 0 to 5, indicating how many RP thresholds were exceeded. A Cox regression approach was used to calculate the hazard ratio associated with the number of profiles above threshold. A final Kaplan–Meier analysis was performed between patients scoring 0–2 and patients scoring 3–5.
The imaging data from the final time point before death were processed as described above in 2 additional patients undergoing autopsy. These patients underwent treatment. The patients' brains were sectioned in the same orientation and thickness as their last MRI. Tissue specimens were acquired from areas suspicious of tumor, including regions highlighted by RPs, then paraffin-embedded and stained using hematoxylin and eosin. The histology was then coregistered to the T1 image using methods previously established (24, 28). To assess the underlying pathology, the RP was overlaid on the histology.
Five profiles were identified as highly correlated with survival independent of the tumor volume, where higher volumes of each RP were associated with poorer prognosis. These included 2133 within FTU ROIs and 1133, 1213, 1233, and 3133 within T1E ROIs (P < .0005). Figure 2 shows the profiles significantly correlated with survival overlaid on FLAIR (2133 in FTU) or T1+C (all others).
Volume thresholds associated with each RP are shown in Table 1. The number of RPs above threshold predicted patient survival, where patients above threshold in 3–5 of the RPs showed decreased OS compared with those with 2 or fewer RPs above threshold (P < .001), and survival curves associated with prognostic score and risk group can be seen in Figure 3. The addition of each RP above volume threshold was associated with a hazard ratio of 1.44 (P < .001). Profile 1133 was histologically confirmed to contain dense hypercellularity and necrosis in 2 patients. Figure 4 shows the histological validation of profile 1133. Panel A shows a FLAIR image at the time point nearest to death, with profile 1133 mapped onto the patient's image. An outline of the tissue sample is shown overlaid in white. A large cluster of 1133 is within the tissue sample, and panel B shows a magnified view of the area in the tissue sample indicated by the RP.
|Profile||Condition||Cutoff (mm3)||N Above/Below Threshold|
This study presents a method of combining information from a variety of MRIs to create profiles that quantify heterogeneity in tumor appearance. We found that 5 RPs significantly correlated with survival when thresholded by volume. When combined, these profiles created a prognostic indicator that may be used to evaluate prognosis noninvasively before treatment. Figure 3 shows 2 Kaplan–Meier plots, and it represents the risk-stratification model these profiles generate. In the left panel, all 6 prognostic scores have individual Kaplan–Meier curves. Although the sample size is smaller, the survival times stratify on the basis of the prognostic score. The right panel shows a combined prognostic score, where all patients with 0–2 score are grouped together in a low-risk group, and patients scoring 3–5 are considered high risk. There is a statistically significant difference in survival time between low- and high-risk patients (P < .001). At autopsy, hypercellular tumor and necrosis was found in voxels indicated by one of the RPs in 2 patients.
Previous studies have found that the visual appearance of brain tumors on imaging such as mass effect, brain tumor contrast enhancement, degree of T2 edema, and contrast to T2 ratio is related to the genetic phenotype of brain tumors (19). Previous radiomic analysis of MP-MRI has found that imaging features such as standard deviation of energy and gray level run emphasis are significantly correlated with patient prognosis (29). Radiomic analysis typically produces a larger number of agnostic features, features that are mathematically extracted descriptors of tumor properties (14). Image features used in radiomic analysis have been shown to capture distinct phenotypic differences in cancers such as those of the lungs (30–33) and the head and neck (9, 33–36).
One of the goals of this study was to bridge the gap between what radiologists interpret with MP-MRI and the field of radiomics. Agnostic radiomic features such as image entropy or image energy are often found to be associated with a clinical metric of interest. Because these features are typically calculated within an ROI, radiologists are limited to interpreting a score rather than a heterogeneous image. The radiomic profiling method presented here results in a simple-to-interpret 4-digit code that maps tumor heterogeneity in a visually meaningful way. Figure 5 compares 2 patients, one with a small tumor and a high-risk score and the other with a larger tumor and a lower prognostic score. The raw FLAIR and T1+C images, as well as the combined ROI, are shown. The bottom panel groups all profiles into one of 3 categories such that the map is easier to read. Blue indicates profiles not significantly correlated with survival, green represents profiles correlated with survival below threshold, and yellow indicates profiles correlated with survival that are above threshold. Despite the smaller tumor in the patient on the right, 5 profiles above threshold result in a much smaller survival time. Additional analysis of the RP volumetrics provides risk stratification, predicting patient prognosis in a manner that is intuitively understandable to a radiologist. The resulting profile maps may be useful for defining resection margins or targeted radiotherapy.
Current practice in GBM imaging involves regular imaging sessions. Radiomic profiling in its current state best functions at the first scans after diagnosis, potentially directing, radiotherapy or surgical resection. Future studies will validate this method over time, accounting for treatment effects and monitoring if changes in RPs continue to correlate with prognosis. The current method is semiautomatic, involving human input only in checking the ROIs after they are initially identified by thresholding. With a more advanced ROI-detection method, this approach could become a fully automated tool at the physicians' disposal. No additional scans would be required to implement this method. The time to generate the RP maps is under 30 seconds assuming the ROIs are already available.
The 4 MRI contrasts included in this analysis are standard in most clinical brain tumor protocols. ADC, which measures the diffusion of water within tissue, has been shown to be inversely correlated with tissue cellularity (37). An increase in ADC has also been hypothesized to indicate cell death (38). FLAIR is a heavily T2-weighted fluid inversion sequence, which suppresses cerebrospinal fluid, so that excess fluid within the brain becomes hyperintense (39). FLAIR hyperintensity is generally interpreted as a combination of infiltrative tumor cells and vasogenic edema.
The code 1133 was hypothesized to be the most significant profile before analysis because it best describes the appearance of invasive tumor on all the following 4 contrasts evaluated: dark on T1, diffusion restricted on ADC, enhancing on T1+C, and within the region of edema on FLAIR. The environments 1113 and 1123 (the same profile with different T1+C values) are not significant in any condition, suggesting that contrast enhancement was the important characteristic in this profile. Four of the 5 significant profiles take the form XX33, indicating that they are contrast-enhancing and hyperintense in the FLAIR. Interestingly, no profiles outside of the contrast-enhancing lesion were found to be indicative of OS.
The volume threshold for each profile was unexpectedly variable, where some profiles, such as 1133, have a near-zero threshold, whereas others, such as 1233, have a threshold over 450 mm3. Future studies should assess phenotypic differences in these RPs.
There were several sources of potential error in this study. Segmentations of both the T1E and FLAIR hyperintensity were performed semiautomatically, with manual correction of misclassified voxels, and thus, are prone to potential human error. Variable imaging parameters, magnet strength, and vendor may also contribute as sources of error. Because of the retrospective nature of the study, images were acquired over the course of several years under similar but slightly varying imaging parameters. The data set would likely produce more robust profiles if the parameters and imaging systems had been identical. We controlled for this by intensity-normalizing each image.
These profiles will have greater clinical importance if they can be linked to a distinct histological pattern. Each profile may describe a specific phenotype or microenvironment. Further histological validation with more patients is necessary. Treatment effects, however, may change the profile volume, threshold, or which profiles are most useful in predicting prognosis. Monitoring profile growth over time will allow us to better control for both treatment effects and time dependence. A variation in prognostic score over time may provide useful clinical information.
In conclusion, this study presents an easily interpretable method for creating radiographic profiles by combining intensity information from multiple MRI scans. When thresholded by volume, 5 RPs significantly correlated with survival. A prognostic indicator that combined the 5 may potentially be used to evaluate prognosis noninvasively before treatment. We also pathologically validated that voxels indicated by one of the RPs contained hypercellular tumor and necrosis in 2 patients. The method presented in this paper may prove clinically useful in providing a risk-stratification model for clinicians treating newly diagnosed GBM.