Gliomas are an aggressive and often deadly form of primary intracranial tumor, representing 81% of malignant brain tumors (1). Current standard treatment includes surgical removal of the tumor followed by administration of radiation and chemotherapy (2, 3). Critical to maximizing the efficacy of these treatments is detecting interval tumor progression and monitoring areas of potential treatment effect using noninvasive imaging. Multiparametric magnetic resonance imaging (MP-MRI) is central to diagnosing a glioma, monitoring its progression, and assessing the efficacy of treatment, as well as providing information regarding neuroanatomy and the physical properties of the tissue. Typical clinical protocols include pre- and postcontrast T1-weighted images (T1, T1+ C), and T2-weighted fluid attenuation inversion recovery (FLAIR) images, which are used to delineate enhancing, nonenhancing, and necrotic tumor components (4, 5). Diffusion-weighted images are also often collected to calculate the apparent diffusion coefficient (ADC) map, which is used prognostically to identify areas of restricted diffusion that may manifest as tumor recurrence (6–9) or stroke (10–12). Despite the current utility of MP-MRI, the pathological heterogeneity of high-grade brain cancers often confounds traditional imaging signatures, while pathological validation studies limited to using biopsy cores as the ground truth may fail to capture tumor pathology in noncontrast enhancing regions (13). Thus, advancements in clinical imaging and tumor monitoring are necessary to more accurately identify the active, nonenhancing tumor components and inform improved treatment directions.
Recently, MP-MRI has been used for a quantitative feature extraction process known as radiomics. These radiomic features attempt to quantify aspects of an image such as intensity distributions, spatial relationships, and textural heterogeneity (14–17). Feature extraction is typically calculated over a given region of interest (ROI), with shape-based features, histogram-based first-order features, and a range of texture features across transformations of the original ROI matrix (18). Extensions of these feature extraction techniques have calculated similar features for 3-D wavelet decompositions of the original image, which split the original image into its high- and low-frequency components to provide enriched edge-related information (19, 20). The end result of these image processing techniques includes mineable data sets of quantitative information that can be compared with clinically relevant metrics. In studies of gliomas and other brain cancers, several groups have shown the utility of radiomic features in predicting prognostic factors such as survival time (21–23), antiangiogenic treatment response (24), and various molecular and genetic subtypes of gliomas (25, 26).
Despite the promising results shown with standard radiomics-based techniques, several methodological limitations have prevented their widespread adoption. Questions regarding the generalizability of radiomics-based predictive models have been raised by studies assessing the stability of the features across MR acquisition parameters (27–30). Current studies of radiomic features tend to focus on feature extractions across the whole tumor for subject-level classification, but studies of localized radiomic features are few. Both the clinical need for better localization of active, nonenhancing tumor and the desire for biological validation of radiomics-based signatures further highlight this gap in the current literature. An initial step in bridging this gap is to show localized associations between radiomic features and similar features of the underlying histology to demonstrate congruence between MR-derived features and features of the underlying tissue.
To address this initial step toward localized radiomics-based modeling, this study compares localized radiomic features to similarly calculated features of coregistered histology images, referred to here as “histomic” features. This study uses MRI coregistered tissue samples from patients with brain cancer that were acquired at autopsy to explore the association between tile-based radiomic features and their histomic analogs. Specifically, this study tests the hypotheses that radiomic features show substantial associative relationships with their histomic analogs and the relationship between radiomic and histomic features is stable across acquisition field strengths.
Sixteen patients with pathologically confirmed primary brain tumors were enrolled for the brain tissue component of this IRB-approved study. A brief outline of the clinical history of each subject is presented in Table 1, and a diagrammatic representation of the data collection process is presented in Figure 1.
i] Overall survival time is calculated from first surgery, except in the case denoted with *, in which case surgery was not performed and survival is calculated from first appearance on MRI.
ii] Abbreviations: RT, radiation treatmen; TMZ, temozolomide; adj. TMZ, adjuvant temozolomide; Bev, bevacizumab; CCNU, lomustine; VCR, vincristine; Thio-TEPA, thio triethylene thiophosphoramide; INF, interferon; PLDR, pulsed low-dose rate radiotherapy; 13-CRA, isotretinoin; CPT-11, irinotecan; TTF, tumor treating fields.
Image Acquisition and Preprocessing
T1, T1+C, FLAIR, and DWI images were acquired from each subject. For each subject, these scans were the last MRI acquired for clinical purposes before death. Images were acquired on our institutional MRI scanners, including 1.5 T and 3 T GE (General Electric Health, Waukesha, Wisconsin) and 1.5 T Siemens (Siemens Healthineers, Erlangen, Germany) magnets. Acquisition parameters for an example subject at 1.5 T include (repetition time/echo time) the following: T1 spin-echo sequence (T1), 666/14 milliseconds; postcontrast T1 spin-echo acquired with intravenous gadolinium (T1+C), 666/14 milliseconds; ADC), calculated from DWI images acquired with an spin echo-echo planar sequence using b = 0 s/mm2 and b = 1000 s/mm2, 10 000/90.7 milliseconds; and FLAIR, acquired with an inversion recovery sequence, 10 000/151.8 milliseconds, and TI of 2200 milliseconds. All T1-/T2-weighted images acquired with 0.43 × 0.43 mm in-plane resolution (matrix range = 512 × 512 voxels) and a slice thickness of 5 mm, with ADC acquired at 0.86 × 0.86 mm in-plane resolution and resampled to the same resolution as the T1-/T2-weighted images. Two Siemens data sets acquired at 1.5 T and 14 GE data sets (5 acquired at 1.5 T, 9 acquired at 3 T) were included in this study.
All magnetic resonance (MR) images for a given patient were rigidly aligned with the FLAIR image using FSL's FLIRT tool (31–33). After coregistration, T1, T1+C, and FLAIR images were intensity-normalized by dividing voxel intensity by its standard deviation across the whole brain (34). This allowed for comparable scaling across subjects for radiomic features calculated from the qualitative MR imaging. Because ADC is considered a quantitative measurement, values from ADC images were not intensity-normalized. Images were visually inspected before and after processing to ensure images were free of any artifacts that may confound analyses.
Ex Vivo Histology Processing
Upon death, the brains of each subject were removed during autopsy and placed into 10% buffered formalin within a 3D-printed cage based on the subject's most recent MRI to maintain structural integrity with respect to the imaging data during fixation, as previously published (35). After ∼2 weeks, the brain samples were sliced using slicing molds printed to delineate the axial slices from the most recent MRI. Tissue samples were collected from each subject based on enhancement on imaging, suspected presence of tumor, or pathologist-determined diagnostic relevance. The samples were then processed, embedded in paraffin, cut, and stained with hematoxylin and eosin (HE). The full slides were photographed at 40× magnification using an Olympus sliding stage microscope. Matlab 2018b (MathWorks Inc., Natick, MA) was then used to process each individual tile, followed by downsampling and compiling of all other tiles from the sample to generate a single image for each slide. Upon examination of these images, samples from each subject showed representative regions of both tumor and nontumor characteristics, allowing assessment of the radiomic–histomic relationship across a range of tissue pathologies.
Previously published in-house software (Hist2MRI Toolbox, written in Matlab) was used to precisely align histology images to the MRI (21, 35–38). Manually defined control points were applied to align each composite histology image with the analogous anatomical features in the FLAIR sequence MRI. Digital photographs were taken at the time of the brain cutting and sample collection to precisely define the location of the histological sample with respect to the MRI. Samples were spatially aligned to the MRI slice that best represented the sample's location by visually inspecting photographs of the brain slices acquired before and after each sample was collected (21, 35–38).
For each grayscale, MR-space histology image, regions of interest (ROI) were manually drawn to designate areas of the image with valid histological information (ie, free from tears, folds). These ROIs were then used to generate tiles for use in radiomic feature extraction using a 10- × 10-voxel sliding frame. Single-voxel strides in each dimension were used to define tiles across the ROIs (n = 54,067), which were used as localized masks for calculation of the radiomic and histomic feature sets.
Pyradiomics v2.1.1 was used to generate the radiomic and histomic feature sets for each tile (18, 20). First-order features (FO; n = 18), gray-level co-occurrence matrix (GLCM; n = 24), gray-level dependence matrix (GLDM; n = 14), gray-level run length matrix (GLRLM, n = 16); gray-level size zone matrix (GLSZM, n = 16), and neighboring gray tone difference matrix (NGTDM, n = 5) features were calculated for each MR image. The same features were additionally calculated on eight 3D wavelet decomposition (3DWD) images of each MRI, generated by applying all combinations of high- and low-pass filters in each dimension. This resulted in a total of 837 radiomic features per MRI modality. The same 837 features were then calculated for the grayscale histology image at MR-resolution for each tile, resulting in an analogous histomic feature set.
Experiment 1: Correlative Analyses.
To examine direct monotonic associations between analogous radiomic–histomic feature pairs, Spearman's rank correlations were computed between the radiomic feature from each MRI and its histomic analog. Peripheral analyses revealed that time between MR acquisition and death influenced the radiomic–histomic relationship (see online supplemental Figure 1); therefore, results are reported both before and after correction for this effect using partial Spearman correlations. Owing to the large tile sample size of this analyses relative to the overall subject count, P-values were not considered a valuable indicator of meaningful effects (eg, a conservative Bonferroni correction for the 3348 comparisons in each heatmap gives an alpha threshold of 0.000015, which corresponds to an indiscriminatory critical value of ρ = 0.0186); thus, results are reported solely in terms of effect size.
Experiment 2: Stability Analyses.
Stability analyses were conducted to observe the effects of acquisition field strength on the generalizability of radiomic–histomic relationships. The tilewise log Euclidean distance (TLED) between the overall radiomic feature set of each MRI and the overall histomic feature set was used to assess global differences in radiomic–histomic associativity. The TLED probability density function for 1.5 T and 3 T GE scans was plotted to assess the global effects of field strength. Differences in distributions were quantified using the Kolmogorov–Smirnoff (K-S) test. In addition, feature-wise assessments of vendor and field strength effects were performed using linear mixed-effect models. Separate models were fit to assess these effects on the normalized difference between each analogous radiomic–histomic pair. The acquisition field strength was included as a main effect in each model, with time between MRI and death included as a covariate as well as nested random effects of subject and tissue sample.
Additional assessments of the effect of MRI scanner vendor within 1.5 T scans were also conducted, mirroring the analysis procedure for acquisition field strength (TLED distribution assessment and mixed-effect models). Although these results provide a quantification for vendor effects in this sample, the very small number of Siemens scans severely reduces the generalizability of these analyses. Therefore, these results have been included as supplementary material (see online supplemental Figure 2).
Experiment 1: Correlative Analyses
Figure 2 shows the resulting correlation heatmaps, displaying the Spearman's ρ between the radiomic and histomic versions of each feature. FO features tended to show the strongest radiomic–histomic association, with very few GLCM, GLDM, GLRLM, GLSZM, and NGTDM features showing meaningful radiomic–histomic relationships. Across the 3DWD images, the strength of association between FO features tended to increase, whereas higher-order associations tended to dissipate across the wavelet decomposition images. Radiomic–histomic associations sorted by strength are presented in Figure 3 to compare feature strengths across image types, with T1+C and FLAIR images showing the strongest radiomic–histomic relationships. Overall, controlling for time between MR acquisition and death tended to decrease the strength of radiomic–histomic associations, but did not affect the general trends regarding image types and feature sets seen in the uncorrected data.
Experiment 2: Stability Analyses
TLED distributions for each acquisition field strength are presented in Figure 4. Most scans displayed minor-to-moderate influence of field strength on the overall TLED, mostly in the form of shifts in centrality. The feature-wise assessment of these confounds is presented in Figure 5, which plots the standardized mixed-model coefficient of the confound against the radiomic–histomic association for each feature. These plots indicate that the features showing the strongest radiomic–histomic relationships tend to have less substantial confounding effects across most image modalities. A statistical summary of the top 10 highest associated features is presented in Table 2.
i] Spearman correlations (ρ) are presented before and after covarying time between MRI and death, and mixed-model coefficients (β) are given for the effect of field strength on the radiomic–histomic relationship. First-order features account for all 10 highest-ranked features presented for each image, particularly those calculated from 3DWD images.
This study explored the histological underpinnings of tile-based MP-MRI radiomic features in patients with GBM and other brain cancers. Several radiomic features showed substantial associations with their histomic analogs (ρ > 0.2), suggesting that these aspects of the MRI directly characterize the same features of the underlying tissue histology. These features are shown to be relatively robust across acquisition field strength, a potential confound to the generalizability of these relationships. These findings, taken as a whole, begin to provide a neuroanatomical context for radiomics-based models of brain cancer characteristics.
To the best of our knowledge, this study is the first assessment of localized relationships between radiomic features and analogous histomic features calculated on histology samples. Our findings suggest that radiomic features are able to capture localized information about the underlying histology of the tissue, motivating a thrust toward radiomics-based models for localized tissue information. Past studies attempting to map localized tissue information using MRI have focused on the use of biopsy cores as the source of ground truth (39–42). Although this allows for characterization of tumor in the MRI-enhancing region, the use of autopsy samples in this study enabled an assessment of the radiomic–histomic relationship in both tumor and nontumor tissue. Whole-brain generalizability will improve the clinical utility of MRI-based tissue feature maps. Thus, this study motivates future investigations using autopsy samples as ground truth to capture whole-brain heterogeneity of tissue features.
Despite the benefits of measuring the radiomic–histomic relationship across autopsy samples, the time between the MRI scan and tissue fixation at death becomes a new confounding factor, as subtle changes in the disease state may occur during this time. When assessing radiomic–histomic relationships after controlling for the duration of this period, nearly all features show mild-to-moderate decreases in associative strength. This decline suggests that modest correlations between radiomic/histomic features and time between MRI and death may gently inflate the associations seen in the uncorrected data, likely due to within-subject variance given the small sample size studied here. Despite these effects, the overall strength and patterns of association remained largely intact, suggesting that minor statistical artifacts in the uncorrected data did not falsely manifest these general trends. Larger studies are warranted to better assess the validity of this claim, as well as to provide better characterizations of MRI-to-fixation duration effects on different MRI-to-tissue mappings.
The most notable contrasts within the associative analyses (Figure 2) include the general strength of FO features compared with that of various higher-order features, as well as the modulation of this relationship across the 3DWT. Generally, FO features presented with more substantial associations than those of higher-order features, although higher-order features of the T1+C and FLAIR images still showed some weak associations. This split in associative strength was compounded by the 3DWT results, with FO features increasing in associative strength and higher-order features decreasing in strength across most wavelet decompositions. The lower-frequency wavelet decompositions (LLL and LLH) deviated from this pattern and showed some increase in higher-order associations. Future research into the interfeature relationships will provide insight as to whether these wavelet decompositions provide new information or simply preserve higher-order associations found in the original image.
Across images, FO total energy frequently showed the strongest radiomic–histomic associations. This feature is calculated as the sum of squares across the intensity values of the tile, and then scaled by the volume of the tile. The consistent associations seen for this feature suggest that areas of enhanced contrast on the MRI associate with intensity changes in the histology, which could reflect differences in tissue composition and presence of tumor pathology. Other FO features with substantial associations include measures of first-order spread such as range, variance, and mean absolute deviation. These features generally characterize the degree of intensity variance across each tile, which may suggest that MRI features are able to capture edge-related information in the histology. The enhanced first-order associations seen across the edge-enhancing 3DWT support this notion as well, generally suggesting that first-order intensity and spread characteristics of the histology are most strongly preserved across MP-MRI.
Comparing across scan types, the FLAIR image tended to have the strongest associations with histology texture (Figure 3). Given the immense diagnostic utility of FLAIR images in delineating different tumor regions, this does not come as a great surprise, and it supports our use of the FLAIR image as the optimal reference image for histology coregistration. The T1+C image also showed substantial radiomic–histomic associations, which suggests that the addition of contrast agent enhances the degree of texture preservation seen in the MRI. The features of the T1+C image with the strongest radiomic–histomic relationships also tended to be more stable across acquisition field strength than those of the FLAIR (FLAIR TLED D = 0.447 vs T1+C TLED D = 0.089). These results suggest that the contrast agent may provide new information relevant to tissue texture beyond that of the standard T1 image, which showed substantially weaker radiomic–histomic associations.
The weak radiomic–histomic associations observed in the T1 and ADC images relative to the FLAIR and T1+C images (Figure 3; lower feature associations across nearly all ranks) imply that these MRI modalities are not as representative of underlying tissue histology texture. This does not necessarily discount the predictive utility of these images as much as it highlights what sort of information these modalities may provide. Whereas the magnitude of the ADC image in a given region may provide robust and insightful information regarding some diagnostic criteria, the features of the ADC and T1 images may not be as relevant to providing histological texture information. This distinction is important as studies move toward assessing the biological underpinnings of radiomics. The lack of texture preservation seen in these results points future studies toward nontextural candidates for biological validation of radiomic features derived from these modalities. Given that ADC is thought of as a quantitative metric, the influence of field strength on some radiomic–histomic relationships was unexpected. This is likely due to subject-specific differences that manifest as group differences in our small sample; however, this could also reflect an amplification of subtle field strength differences in ADC values owing to tile-based sampling and higher-order feature calculations. Further, larger-scale research into how radiomic features behave using different acquisition parameters, field strengths, and scanner vendors will likely provide insight toward the mechanism for this effect.
Our clinical imaging data sets contained a small range of acquisition field strengths, which allowed for preliminary assessments of how these factors influence the radiomic–histomic relationship. Previous studies have shown several radiomic features susceptible to differences across different scanners (43–45); this study adds to this literature by providing a breakdown of how well features represent histological texture across different acquisition field strengths. As desired, the features with the greatest degree of association between radiomic and histomic analogs tended to have the least substantial effects of field strength (Figure 5, lower β values for higher ρ values), with the notable exception of the evenly confounded FLAIR features (Figure 5, even spread of βs across ρ values). Statistical artifacts of the confounds may influence the weaker associations observed in this data set, but generally these results point to a subset of features characterized by both strong and stable radiomic–histomic associations. Larger follow-up studies assessing radiomic–histomic relationships may be able to reveal additional features that accurately characterize tissue texture for particular scanner vendors and acquisition parameters.
This study is not without its limitations, however. The relatively small sample size of this study calls for replication of these results in larger data sets to confirm the generalizability of the patterns seen here. This study was also limited to patients with primary brain cancer, which reduces the scope of the radiomic–histomic relationship these results imply. Future studies of radiomic–histomic relationships in other neurological diseases will be an important step in establishing more general cases for the relationship between MRI and tissue texture. The use of a tile-based prediction to study localized information does not allow for the use of shape- and size-based radiomic features, which often provide useful weights in radiomics-based modeling of prognostic factors. The tile-based extraction method also complicates quantitative feature mapping owing to uneven sampling of center and edge voxels. Future studies studying voxel-level features may be able to provide these quantitative maps to visualize association patterns across different brain regions, as well as across different tumor components.
Although this study statistically controlled for the duration between MRI acquisition and death, it is possible that these effects are not adequately addressed as covariates in this small sample. Future, large-scale studies of autopsy data will be essential to characterize the magnitude of this effect, as well as to address optimal strategies to account for these effects. This study focused on assessing stable feature relationships across all subjects, but specific radiomic signatures of treatment effects and treatment resistance may also be present in ways this treatment-heterogeneous sample could not reveal. Prospective imaging studies may aid in assessing radiomic–histomic relationships across different periods of treatment, and thus may be better suited to study these imaging signatures. Animal studies may also be able to provide insight into the radiomic–histomic relationship by allowing for control over time between imaging and tissue resection, as well as for larger-scale studies.
In conclusion, this study presents a novel characterization of the radiomic–histomic relationship in an attempt to provide a neuroanatomical basis for radiomics-based analyses. These results show a substantial degree of heterogeneity in the strength and stability of radiomic–histomic relationships but reveal a subset of radiomic features that stably reflect information about the underlying histology. This study provides the groundwork for future investigations into the quantitative pathological validation of radiomics analyses as currently performed, as well as underscores an opportunity for radiomics-based predictions of localized histological features with diagnostic and prognostic utility.