Abstract
Background This study identified a clinically significant subset of glioma patients with tumor outside of contrast-enhancement present at autopsy, and subsequently developed a method for detecting non-enhancing tumor using radio-pathomic mapping. We tested the hypothesis that autopsy-based radio-pathomic tumor probability maps would be able to non-invasively identify areas of infiltrative tumor beyond traditional imaging signatures.
Methods A total of 159 tissue samples from 65 subjects were aligned to MRI acquired nearest to death for this study. Demographic and survival characteristics for patients with and without tumor beyond the contrast-enhancing margin were computed. An ensemble algorithm was used to predict pixelwise tumor presence from pathological annotations using segmented cellularity (Cell), extracellular fluid (ECF), and cytoplasm (Cyt) density as input (6 train/3 test subjects). A second level of ensemble algorithms were used to predict voxel-wise Cell, ECF, and Cyt on the full dataset (43 train/22 test subjects) using 5-by-5 voxel tiles from T1, T1+C, FLAIR, and ADC as input. The models were then combined to generate non-invasive whole brain maps of tumor probability.
Results Tumor outside of contrast was identified in 41.5 percent of patients, who showed worse survival outcomes (HR=3.90, p<0.001). Tumor probability maps reliably tracked non-enhancing tumor in the test set, external data collected pre-surgery, and longitudinal data to identify treatment-related changes and anticipate recurrence.
Conclusions This study developed a multi-stage model for mapping gliomas using autopsy tissue samples as ground truth, which was able to identify regions of tumor beyond traditional imaging signatures.
1. Introduction
Glial brain tumors are the most common form of central nervous system originating tumors, with a yearly incidence of around 6 per 100,000 persons1. High-grade gliomas such as glioblastomas (GBMs) are often difficult to treat due to their aggressiveness and pathological heterogeneity, contributing to one- and five-year survival rates of 41 percent and 5 percent, respectively2,3. Current standard of care requires precise tumor localization to maximize the efficacy of frontline treatments such as surgical resection and targeted radiation therapy. Studies examining extent of resection in glioma patients have demonstrated maximizing tumor removal improves patient survival outcomes, underscoring the importance of identifying the full extent of gliomas4–7.
Magnetic resonance imaging (MRI) is currently the primary method for non-invasive disease monitoring of gliomas. Gadolinium contrast agents are used to highlight disruptions in the blood brain barrier due to tumor angiogenesis, which results in an enhanced signal in T1-weighted imaging that is used to define the primary tumor mass8–10. Hyperintensities on T2-weighted fluid attenuated inversion recovery (FLAIR) images correspond to a mixture of infiltrative tumor and edema, though differentiating between tumor and non-tumor remains a source of uncertainty in imaging interpretation11–14. Apparent diffusion coefficient (ADC) images derived from diffusion-weighted imaging have shown promise in highlighting areas of diffusion restricting hypercellularity associated with tumor, though recent studies validating this signature beyond the contrast-enhancing margin have disputed the strength of this relationship14–19. Machine learning approaches have also sought to maximize the amount of clinically relevant information we can extract from these non-invasive images, employing recent advances in computing to automatically segment the radiologist-defined margin, identify patient-level genetic signatures such as IDH1 mutation and MGMT methylation status, and predict cellular-level information using biopsy samples as validation20–24.
Despite the promise of these recent techniques, identifying areas of tumor beyond the contrast- and FLAIR-enhancing margins remain particularly difficult. Studies of autopsy tissue samples aligned to MRI acquired near death have found areas of active tumor as far as 10 cm beyond the treated margin, indicating the need for improvements in tumor tracking17,25,26. However, identifying regions of tumor beyond current signatures requires ground truth tissue samples from beyond the suspected tumor region, which is beyond the capabilities of studies using biopsy tissue and radiologist annotations as the gold standard. Therefore, studies with access to tissue beyond contrast enhancement are essential for identifying the true extent of tumor, particularly in the less-studied post-treatment state. This study uses autopsy tissue samples, collected both within and beyond the traditional tumor margin, to assess the prevalence of tumor outside the MRI-defined margin and develop a multi-stage model for non-invasive tumor probability mapping. Specifically, we tested the hypothesis that an MRI-based model for tumor probability trained on autopsy tissue samples can track glioma invasion beyond the currently defined margin.
2. Methods
2.1. Patient Population
This study was approved by the Institutional Review Board of Medical College of Wisconsin. Written, informed consent was obtained from 65 patients to participate in this study, each diagnosed with a glioma in concordance with the 2021 WHO classification standards for brain tumors. These participants were divided into a training dataset used for model development and a test set used to test model generalizability. Clinical and demographic characteristics of this sample, as well as for the subsamples relevant to this study, are presented in Supplemental Table 1. A schematic of the overall study outline and data processing steps is presented in Figure 1.
2.2. MR Image Acquisition and Preprocessing
Clinical imaging from each patient’s last MRI session prior to death was collected for use in this study. Due to the use of clinical imaging, acquisition parameters were not standardized across patients. T1, T1+C, FLAIR, and ADC images were selected as the primary focus for this study, as they provide a wide range of information regarding tissue characteristics while remaining ubiquitous across clinical acquisitions. T1, T1+C, and ADC images were each rigidly aligned to the FLAIR image using SPM12 (https://www.fil.ion.ucl.ac.uk/spm). Qualitative images (T1, T1+C, and FLAIR) were then divided by the intensity standard deviation within the brain to normalize values across patients8,14. For the test set, MRI sessions from the patient’s entire clinical history were collected, aligned to the FLAIR nearest to death, and normalized to assess the ability of the model to track tumor growth patterns over time.
2.3. Tissue Segmentation and Processing
A total of 159 tissue samples were collected from participants at autopsy using previously described methods and placed in 66 mm by 49 mm large format cassettes14,17,26. Following hematoxylin and eosin (HE) staining, digitization at 40X magnification, and downsizing by a factor of 10, a color deconvolution algorithm was used to separate the images into color channels that correspond to stain optical density (i.e., positive hematoxylin staining, positive eosin staining, background)27. This image was then filtered and processed to extract nuclei, cytoplasm (Cyt,), and extracellular fluid (ECF). Both segmentations were then upsampled to approximate MRI resolution, using 50×50 superpixels to convert segmentation masks to cell density (cells/mm2) and ECF/Cyt proportion (ranging from 0 to 1). Additionally, 33 tissue samples from a subset of 9 participants were annotated for presence of infiltrative tumor (IT), tumor with pseudopalisading necrosis (PN), non-tumor necrosis (NtN), and unlabeled tissue (UL) by a pathologist-trained technician.
2.4. MRI-Histology Co-registration
Tissue samples were aligned to each participant’s FLAIR image using previously published in-house software, written in Matlab14,17,26,28,29. Manually defined control points were used to identify architectural landmarks on both the tissue data and the MRI using photos from autopsy to guide accurate placement. These control points were then used to compute a non-linear transform that warped tissue to the same space as the MRI, and regions of interest were drawn to exclude areas of tissue distortion (i.e., rips, tears, folds) and MRI artifacts. Following MRI-tissue co-registration, all MRI and tissue data was resampled to 0.4397 × 0.4397 mm per pixel resolution to harmonize all subjects to the most common acquisition parameters.
2.5. Assessing prevalence of tumor beyond MRI-defined margin
Aligned autopsy tissue was then compared to the T1+C MRI to estimate the prevalence of tumor outside of contrast enhancement (TOC). Tissue samples were visually inspected for tumor presence using histological characteristics (i.e. hypercellularity, mitotic events, pseudopalisading necrosis), as well as any annotations or immunohistochemical staining available (i.e., Ki-67, CD31). Patients were then coded for presence of TOC by comparing manually segmented maps of contrast to the aligned tissue samples (see Supplemental Figure 1 for examples). Chi-squared tests were also used to determine differences in TOC frequency amongst patients who did and did not receive radiation with temozolomide (Rad+TMZ), Bevacizumab (Bev), and tumor treating fields (TTFields), as well as amongst diagnostic grade groups. A Cox proportional hazards regression was then fit to compare survival durations between patients with and without TOC, including time between MRI and death as a covariate. Survival curves were limited to patients who had received standard treatment (Rad+TMZ) to exclude a small number of untreated patients who tended to survive shorter than similar patients who had received treatment. Analysis was performed for the full data set (N=65), as well as on subjects with scans less than 90 days away from death (N=47) and patients with a primary GBM at first surgery (N=37) to compare results across potential confounds.
2.6. Predicting tumor presence from pathological segmentations
A chart indicating the flow of patients across the multiple model training stages is presented in Supplemental Figure 2. The first component of the multi-stage tumor prediction model involved predicting tumor annotations from the segmented pathology data. Several different candidate models were assessed, including k-nearest neighbors, naïve bayes, decision trees, and random-under-sampling-boosted random forest (RUS Tree) models. Each model was trained using pixelwise cell, Cyt, and ECF density as input to predict tumor annotation (IT and PN) versus non-tumor (NtN and UL). Models were trained on slides from 6 subjects and validated on the 3 remaining subjects to assess generalizability. Model performance was evaluated using receiver operator characteristic (ROC) plots and area under the curve (AUC) metrics. The model with the highest AUC was then incorporated into the multi-stage model.
2.7. Predicting pathological segmentations from MRI data
The second stage of the multi-stage model focused on training separate models to predict cell, Cyt, and ECF density using MRI data. Bootstrap aggregating (bagging) random forest models were trained to predict voxel-wise densities using 5 by 5 voxel tiles from the T1, T1+C, FLAIR and ADC as input. This model framework was selected based on our previous publication, which developed a proof-of-concept model for predicting cellularity using MRI data in a smaller sample of patients25. These models were developed on 2/3rds of the full dataset (N=43) and tested on the held-out test set (n=22) to assess generalizability. None of the patients in the pathological tumor prediction test set overlapped with the training set for these radio-pathomic models, ensuring that the test set contained only data unseen in training for every component model. Quantitative performance was evaluated using root-mean-squared error (RMSE) estimates within each test set subject, and example segmentations were plotted and compared to ground truth segmentations to ensure accurate spatial localization of cell, ECF, and Cyt trends. Predicted maps for cellularity were set to range from 0 to 6000+ cells/mm2, which offered good contrast between hypercellular tumor and normal tissue when examining tissue samples, and ECF/Cyt density maps were set to range from 0 to 1, reflecting the minimum and maximum possible values of the measurement.
2.8. Predicting tumor presence from MRI data
Following model training, the best performing pathological tumor prediction model was used to convert whole brain cell, ECF, and Cyt density maps to tumor probability maps (TPMs). These maps were then plotted against predicted segmentations to assess patterns in what the model classifies as tumor, as well as against ground truth tissue segmentations and annotations to assess accuracy at identifying novel tumor beyond traditional signatures. TPMs were then applied to longitudinal images from all applicable timepoints across each test set patient’s clinical history to assess the ability for our tumor prediction model to identify areas of pre-angiogenic recurrent tumor prior to contrast enhancement. Additionally, whole brain TPMs were applied to two external data sets to assess the quality of the maps (i.e., free from artifacts, noise) and ability to identify non-enhancing tumor on MRI data acquired entirely independently of this study. The first external subject was a 61 year old male diagnosed with a recurrent GBM, whose scan was collected by our external collaborators at the University of California – Los Angeles. The second external case was a pre-treatment acquisition from a 50 year old male diagnosed with GBM, taken from the publicly available TCGA-GBM dataset (https://cancergenome.nih.gov/)30,31.
3. Results
3.1. Prevalence of non enhancing tumor
Figure 2 shows the prevalence estimates for TOC presence, as well as corresponding Kaplan-Meier plots for survival analyses. TOC was observed in 41.5 percent of patients, with increased presence amongst GBM (GIV) patients (χ2=10.73, p=0.005), patients who have been treated with Rad+TMZ (χ2=3.99, p=0.046), and patients treated with bevacizumab (Bev) (χ2=5.41, p=0.020). Treated patients with TOC showed decreased survival rates compared to patients without TOC (HR=3.90, 95 percent CI = 2.50-5.30, p<0.001). Patients with MRI data from within 90 days of death also showed increased TOC presence amongst Rad-TMZ and bevacizumab-treated (χ2=4.50, p=0.033 and χ2=4.50, p=0.033, respectively) patients, and reduced survival associated with TOC (HR = 3.85, 95 percent CI = 2.33-5.37, p = 0.001). Primary GBM patients did not show significant differences in TOC presence with respect to treatment but did show reduced survival within patients with TOC (HR = 2.95, 95 percent CI = 1.32-4.57, p = 0.025).
3.2. Predicting tumor presence from pathological segmentations
Figure 3 shows the performance results for the pathological tumor prediction model, along with example tumor probability maps and corresponding ground truth annotations. The RUS Tree model gave the best ROC AUC of 0.857 and was therefore selected as the model used to convert MRI-based pathological segmentations into tumor probability maps. Example predictions show good concordance between areas of high tumor probability and actual tumor presence within both IT and PN tumor types, while correctly avoiding edema and other non-tumor areas.
3.3. Predicting pathological segmentations from MRI data
Segmentation prediction results are shown in Figure 4, along with example whole brain predictions on test set subjects. Each radio-pathomic model had a mean subject-level RMSE value within a standard deviation of the ground truth (cell density Std. RMSE = 0.756, Cyt Std. RMSE = 0.917, ECF Std. RMSE = 0.941), indicating satisfactory model performance for each tissue type. Example predictions indicate regions where these maps improve pathological interpretations of the MRI data, such as discrimination between areas of edema and hypercellularity within the FLAIR hyperintense region and areas of hypercellularity beyond contrast enhancement.
3.4. Predicting tumor presence from MRI data
An example TPM generated from the full multi-stage model are presented in Figure 5. Test set TPM with pathological annotations available showed good correspondence between areas of predicted and actual tumor, with separation between areas of PN (high cellularity, high ECF) and areas of non-PN tumor (high cellularity, normal ECF) seen on the corresponding segmentation maps. Tumor was also accurately observed beyond the contrast enhancing region. External data and longitudinal TPMs are presented in Figure 6. Predictions on external data sets demonstrated that TPMs generated from this framework provided low-noise, interpretable maps on new data, as well as identified new areas of possible tumor infiltration beyond contrast enhancement. Longitudinal predictions for the example subject showed the TPM in response to treatment, where areas of heightened tumor probability tended to fade after administration of Bev+TMZ and returned following treatment cessation. These maps also identified tumor crossing-over into the contralateral hemisphere as early as 374 days prior to death despite a lack of contralateral contrast enhancement across the patient’s imaging. This area of crossing over was confirmed to be tumor at autopsy.
4. Discussion
This study developed a multi-stage tumor prediction model that accurately identifies areas of infiltrative tumor beyond traditional imaging signatures. A pathological tumor prediction model based on histology segmentations was combined with MRI-based maps of tissue composition to bridge the gap between histopathological assessment and non-invasive tracking of glioma. We validated this technique using autopsy tissue samples collected near death as ground truth, and successfully applied our models to external data to demonstrate generalizability. Furthermore, we showed that these models can detect areas of tumor recurrence prior to contrast enhancement in retrospective longitudinal data from our test set patients, indicating further clinical utility in disease tracking across patients’ clinical histories. This study is the first of its kind to provide tumor probability maps validated using tissue well-beyond the current treatment margin, and is, to date, the largest study of radio-pathomic signatures at autopsy in brain cancer.
This study demonstrates that around half of gliomas exist in part beyond the treated margin, with higher rates observed in higher grade tumors such as glioblastomas. Since this study only used select tissue samples from areas of suspected tumor that could be aligned accurately to the MRI after tissue processing, this frequency is likely a lower bound, as it is possible that areas of tumor are missed in the tissue sample selection process. However, this prevalence estimate highlights that tumor outside of contrast enhancement is present in most cases, and targeted future MRI-based research should focus on these regions specifically to improve non-invasive tumor localization. In addition, future research at the molecular level could identify pathologic or genetic signatures of enhancing versus non-enhancing tumors, because a better understanding of how tumor evades current detection mechanisms could drive future treatments directed at these non-angiogenic areas of tumor. This is particularly crucial to understanding tumor in the treated state, as our preliminary assessment of treatment effects on TOC/TOF highlight treatment’s impact on radio-pathomic relationships. Bevacizumab in combination with Rad+TMZ treatment was associated with greater tumor outside of contrast compared to bevacizumab-naïve patients, which is hypothesized to be caused by the anti-angiogenic treatment disrupting the co-occurrence of tumor and contrast enhancement. In the future, more targeted research should be conducted to elucidate the mechanism of this relationship further, as well as to understand the influence of factors such as treatment duration and timing on non-enhancing tumor invasion.
To detect regions of tumor beyond traditional imaging signatures, we developed a multistage model to predict tumor using autopsy tissue samples as ground truth. Prior research attempting similar tasks have used biopsy tissue samples, which, while useful, inherently fail to capture areas of tumor beyond the treated margin19,32–34. Our prior studies of autopsy tissue samples have highlighted how imaging signatures validated within contrast enhancement, such as an inverse relationship between ADC and cellularity, tend to show much weaker relationships when validated using tissue from beyond the enhancing region14,17. Our prior paper also developed a proof-of-concept model for predicting cellularity using autopsy tissue samples, which, while insightful, does not directly correlate with areas of tumor, as not all tumor regions (i.e., areas of pseudopalisading necrosis) are strictly hypercellular and not all areas of increased cellularity (i.e., areas of immune response) are tumor25. This study sought to expand upon this method by incorporating multiple tissue segmentations to predict tumor presence itself. In addition to providing precise assessments of tumor presence to distinguish between regions of tumor and non-tumor, the tissue segmentations for cellularity, ECF, and cytoplasm aided in interpreting different features of tumor, such as distinguishing between areas of PN and non-PN tumor, as well as differentiating between areas of tumor and vasogenic edema within FLAIR hyperintensity. This highlights how providing a rich pathological ground truth can vastly improve non-invasive tumor tracking, particularly using machine learning techniques to map relationships invisible to the naked eye.
In addition to applying our model to a held-out test set to assess model performance near death, we applied our model to longitudinal imaging throughout patients’ clinical history to assess our model’s ability to predict recurrence prior to contrast enhancement. A subset of cases showed areas of predicted tumor that became contrast-enhancing tumor months later, and some subjects even showed areas of infiltrative tumor beyond contrast enhancement months prior to death that never became contrast-enhancing but did in fact contain tumor at autopsy. Though future research using biopsy samples is required to confirm tumor presence at timepoints prior to death, these results suggest that TPMs may identify tumor prior to its visual appearance on traditional imaging. By applying this technique to retrospective data and previously collected clinical trials, these maps may be able to reveal new signatures of treatment response in terms of reduced hypercellular volume and low tumor probability, which in turn could identify subsets of patients that selectively respond to treatment.
4.1. Limitations
Due to the use of autopsy tissue samples, the time between imaging and death is an important limiting factor in this study. It is likely that tumor continues to grow near death, and imaging acquired months before death may not reflect how the patient’s imaging would look immediately prior to death. We have accounted for this in our study by restricting our frequency assessment of tumor outside contrast enhancement to imaging within two months of death and have matched our train and test data sets for time between MRI and death. However, future studies using longitudinal imaging to model tumor growth rates may control for this factor more precisely. Post-mortem MRI has also been considered as an alternative, however the changes in tissue properties post-fixation and inability to collect common brain cancer imaging acquisitions (i.e., T1+C) would hinder the ability for our models to predict tumor on clinical data. It is also possible that tissue distortions that have occurred during processing may result in imperfect alignment with MRI data; however, our previously published tissue processing protocol controls for this factor at every step of the tissue collection process, and regions of interest defined after non-linear warping are used to exclude regions of suboptimal tissue quality. Additionally, while this is the largest study of its kind using an unprecedented amount of autopsy tissue data aligned to MRI, the validation set for this study was relatively small compared to other machine learning studies. Our predictions on external data sets highlight the feasibility of applying this method to large data sets for further validation, but precise validation using tissue samples remains an area of future research beyond the scope of this study.
4.2. Conclusion
This study used autopsy tissue samples aligned to MRI to develop a multi-stage model for tumor probability to predict areas of non-enhancing tumor that occur in most glioma patients in our data set. This model accurately identified regions of tumor beyond traditional imaging signatures, provided information for distinguishing between hypercellular and necrotic areas of tumor, and identified distinct phenotypes based on the maps that differ in survival and treatment efficacy. This technique has the potential to improve disease tracking in the post-treated state and could expand the treated margin to encompass additional infiltrative tumor missed by current tracking techniques. This improvement in disease tracking can thus improve treatment efficacy and clinical outcomes.
Data Availability
The MRI and histology data sets used in this study are not currently publicly available due to data usage agreement restrictions. De-identified data is available upon reasonable request.
5. Disclosures
Authors PSL and JMC receive funding from Novacure Inc. unrelated to the scope of this study.