Abstract
Glioblastoma is the most common malignant primary brain tumor with poor overall survival. Magnetic resonance imaging (MRI) is the main imaging modality for glioblastoma but has inherent shortcomings. The molecular and cellular basis of MR signals is incompletely understood. We established a ground truth-based image analysis platform to coregister MRI and light sheet microscopy (LSM) data to each other and to an anatomic reference atlas for quantification of 20 predefined anatomic subregions. Our pipeline also includes a segmentation and quantification approach for single myeloid cells in entire LSM datasets. This method was applied to three preclinical glioma models in male and female mice (GL261, U87MG, and S24), which exhibit different key features of the human glioma. Multiparametric MR data including T2-weighted sequences, diffusion tensor imaging, T2 and T2* relaxometry were acquired. Following tissue clearing, LSM focused on the analysis of tumor cell density, microvasculature, and innate immune cell infiltration. Correlated analysis revealed differences in quantitative MRI metrics between the tumor-bearing and the contralateral hemisphere. LSM identified tumor subregions that differed in their MRI characteristics, indicating tumor heterogeneity. Interestingly, MRI signatures, defined as unique combinations of different MRI parameters, differed greatly between the models. The direct correlation of MRI and LSM allows an in-depth characterization of preclinical glioma and can be used to decipher the structural, cellular, and, likely, molecular basis of tumoral MRI biomarkers. Our approach may be applied in other preclinical brain tumor or neurologic disease models, and the derived MRI signatures could ultimately inform image interpretation in a clinical setting.
SIGNIFICANCE STATEMENT We established a histologic ground truth-based approach for MR image analyses and tested this method in three preclinical glioma models exhibiting different features of glioblastoma. Coregistration of light sheet microscopy to MRI allowed for an evaluation of quantitative MRI data in histologically distinct tumor subregions. Coregistration to a mouse brain atlas enabled a regional comparison of MRI parameters with a histologically informed interpretation of the results. Our approach is transferable to other preclinical models of brain tumors and further neurologic disorders. The method can be used to decipher the structural, cellular, and molecular basis of MRI signal characteristics. Ultimately, information derived from such analyses could strengthen the neuroradiological evaluation of glioblastoma as they enhance the interpretation of MRI data.
- diffusion tensor imaging
- glioma
- light sheet microscopy
- magnetic resonance imaging
- relaxometry
- tissue clearing
Introduction
Glioblastoma (GBM) is the most common primary brain tumor with poor overall survival (Wen and Kesari, 2008), characterized by neoangiogenesis, invasive tumor growth, and an immunosuppressive tumor microenvironment. MRI is the most important imaging modality for the initial diagnosis and treatment monitoring of a glioma but has inherent limitations in sensitivity and specificity because of a lack of specific contrast agents, MR sequences, and mesoscopic resolution (∼1 mm isotropic resolution at clinical field strength). Therefore, characterizing different areas of malignancy, defining exact tumor borders, and differentiating tumor progression from treatment-related changes such as pseudoprogression remain challenging using established multiparametric MRI techniques. Current research focuses on the potential of innovative advanced MRI techniques such as chemical exchange saturation transfer imaging (Zhou et al., 2022) or MR elastography (Reiss-Zimmermann et al., 2015; Bunevicius et al., 2020) to improve the neuroradiological workup of brain tumors. However, the structural, cellular, and/or molecular foundations of these novel MR signals are incompletely understood.
To overcome these limitations by directly correlating MR changes with ultrastructural findings, we performed MR imaging in healthy animals and three glioma models using high-field preclinical MRI at 9.4 T. Glioma models included syngeneic GL261 gliomas and the patient-derived human glioma lines U87MG and S24 (Osswald et al., 2015) to allow the monitoring of key glioma features including the tumor microenvironment in GL261 tumors, neoangiogenesis in U87MG tumors, and infiltrative tumor growth in S24 gliomas.
To establish a cellular ground truth for the analyses of MRI data, we performed tissue clearing and light sheet microscopy (LSM; “ultramicroscopy”) to obtain three-dimensional datasets of the entire cleared mouse brain. Labeling included intravital dye labeling of the microvasculature using fluorescent lectins (Wälchli et al., 2015; Hahn et al., 2021), tumor cell labeling using fluorescent proteins (S24:tdTomato, U87MG:tdTomato) or myeloid cell labeling in Cx3cr1 reporter mice (Jung et al., 2000). Multiparametric MRI protocols consisted of conventional T1-weighted (T1w) and T2-weighted (T2w) sequences as well as advanced functional sequences such as T2 and T2* relaxometry and diffusion tensor imaging (DTI; Krämer et al., 2019).
MRI and LSM datasets were coregistered to each other and to a reference anatomic atlas (Dorr et al., 2008) of 62 predefined subregions facilitating comparative analyses. Whole tumor and tumor subregions were segmented on LSM datasets, and quantitative MRI parameters were assessed in these regions as well as in atlas-derived anatomic regions. Thus, imaging signatures could be directly related to histopathological LSM features with cellular resolution (axial resolution, ∼5 µm). Our approach can be extended to various preclinical brain tumor or other neurologic disease models. Ultimately, MRI signatures derived from such ground truth-based correlated image analyses could potentially improve patient care by boosting the sensitivity and specificity of MRI data.
Materials and Methods
Animal experiments
All animal procedures were conducted in accordance with the institutional laboratory animal research committee as well as the Society for Laboratory Animal Science guidelines and were approved by the regional administrative authority (Regierungspräsidium, Karlsruhe, Germany; Permit nos. G8/14, G189/12, and G27/17). All animals were housed in cages with four to six mice per cage in a dedicated facility with 12 h light/dark cycle and provided food and water ad libitum. All cell lines used were routinely tested by the multiplex cell contamination test, which detects 28 potential cell culture contaminations (Schmitt and Pawlita, 2009).
Healthy control mice and sham injection
Six 18-week-old healthy female C57BL/6 mice underwent MRI and served as a healthy control group. After that, 2 µl of PBS (Sigma-Aldrich) were stereotactically injected in the right basal ganglia (coordinates: 1 mm rostral and 2 mm lateral from bregma at a depth of 2 mm) using a Hamilton syringe with a fine step motor. Mice were anesthetized with ketamine/xylazine during the procedure. MRI was performed 9 d after PBS injection, and animals served as sham controls to investigate the effects of a needle tract on MRI parameters.
Orthotopic glioma models
GL261.
GL261 is a syngeneic glioma cell line in the immunocompetent C57BL/6J strain (Jacobs et al., 2011), and thus has, in contrast to U87MG- or the S24-models, a normal innate and adaptive immune cell compartment. GL261 glioma cells show a bulky growth pattern. GL261 cells were cultured in DMEM containing 10% FBS, 100 U/ml penicillin and 100 µg/ml streptomycin (all from Sigma-Aldrich). One hundred thousand GL261 cells in 2 µl of PBS (Sigma-Aldrich) were stereotactically implanted in the right basal ganglia (coordinates: 1 mm rostral and 2 mm lateral from bregma at a depth of 2 mm) of five 20- to 30-week-old female transgenic Cx3cr1-enhanced green fluorescent protein (EGFP)/+ C57BL/6 mice that express fluorescently labeled microglia (Jung et al., 2000). Mice were anesthetized with ketamine/xylazine during the procedure. Implantation was performed using a Hamilton syringe with a fine step motor. MRI was performed 18 d after tumor implantation.
U87MG.
U87MG is a human glioma cell line that grows bulky when orthotopically implanted as a xenograft model in immunodeficient mice, and shows intratumoral heterogeneity as well as neoangiogenesis (Breckwoldt et al., 2019). U87MG cells (LGC Standards) were cultured in DMEM containing 10% FBS, 100 U/ml penicillin, and 100 µg/ml streptomycin (all from Sigma-Aldrich). For the in vivo model, 50,000 U87MG:tdTomato cells in 2 µl of PBS were stereotactically implanted in the right hemisphere of five 6- to 8-week-old male athymic NMRI nude mice (Charles River). MRI was performed 28 d after tumor implantation.
S24.
S24 is a patient-derived 1p/19q noncodeleted, IDH (isocitrate dehydrogenase) wild-type GBM stem cell line (Lemke et al., 2012, 2014; Osswald et al., 2015; Venkataramani et al., 2019) exhibiting a diffusely infiltrating growth pattern. S24:tdTomato cells (stably transduced by lentivirus) were cultivated as sphere cultures in DMEM/F-12 under serum-free, “stem-like” conditions including 2% B27 supplement (Thermo Fisher Scientific), 5 µg/ml human insulin (Sigma-Aldrich), 12.8 ng/ml heparin (Sigma-Aldrich), 0.4 ng/ml epidermal growth factor (R&D Systems) and fibroblast growth factor (Thermo Fisher Scientific) as previously described (Venkataramani et al., 2019). Fifty thousand S24 cells diluted in 2 µl of PBS (Sigma-Aldrich) were stereotactically implanted into the right basal ganglia of three 8- to 10-week-old male NMRI nude mice (Charles River). MRI was performed 91 d after tumor inoculation.
MR imaging
MRI was performed on a 9.4 tesla horizontal bore small animal MRI scanner (gradient strength: 675 mT/m; BioSpec 94/20 USR, Bruker) using an 8.4 cm body coil for transmission and a 2 × 2 surface array coil for reception. For imaging, anesthesia was induced with 4% isoflurane (Baxter) in 100% O2 and maintained with 1–1.5% isoflurane in 100% O2 delivered via a nose cone during scanning. The respiration rate was constantly monitored. Animals were placed prone on a Bruker standard MRI bed with an integrated circulating water heating system for maintenance of body temperature. Control mice and mice after PBS sham injection underwent an extensive MR protocol comprising a standard 3D T2w sequence, DTI and T2 and T2* relaxometry. To account for the different characteristics of the three tumor models, the respective MR protocol was adapted to the particular model: in all three models, an identical standard 3D T2w sequence was acquired. For the GL261 cohort, imaging comprised additional DTI and T2 and T2* relaxometry. The protocol of the U87MG cohort included additional T2 and T2* relaxometry, while a gadolinium contrast-enhanced 3D T1w sequence and a DTI sequence were additionally acquired in the S24 cohort. Detailed information on scan parameters is provided in Table 1. Animals of the glioma cohorts were killed after imaging for tissue harvesting and further tissue clearing and LSM analysis.
Detailed MRI sequence parameters
Clearing and light sheet microscopy
Fixation and clearing of mouse brains
For labeling of the microvasculature, animals were injected with fluorescent lectins binding to the N-acetyl-β-d-glucosamine oligomers of endothelial cells (Wälchli et al., 2015; Breckwoldt et al., 2016). Isolectin-FITC or lectin Texas Red dye (12 mg/kg; Sigma-Aldrich or Vector Laboratories) in 100 µl of PBS was injected intravenously 5 min after lectin injection, and animals were killed by a ketamine/xylazine overdose and transcardially perfused with 10 ml of PBS followed by 10 ml of 4% PFA in PBS. Brains were harvested, postfixed with 4% PFA for 24 h, and stored in PBS at 4°C in the dark.
Whole brains of the GL261 cohort were optically cleared using the vDISCO protocol (Cai et al., 2019), which includes nanoboosting steps enhancing the fluorescence signal. In brief, brain samples were incubated in permeabilization solution containing 1.5% goat serum (catalog #16210072, Thermo Fisher Scientific), 0.5% Triton X-100 (Sigma-Aldrich), 0.2% trans-1-acetyl-4-hydroxy-l-proline (catalog #441562, Sigma-Aldrich), 0.5 mm methyl-β-cyclodextrin (catalog #332615, Sigma-Aldrich), and 0.05% sodium azide (catalog #71290, Sigma-Aldrich) in 0.01 m PBS at 37°C for 4–5 d with gentle shaking. Then, the labeling solution was prepared by adding Atto647N-conjugated anti-GFP nanobooster (catalog #gba647n-100, Chromotek) and propidium-iodide (PI; catalog #P4864, Sigma-Aldrich) into the same permeabilization solution at dilutions of 1:600 and 1:250, respectively (8–9 µl of nanobooster and 20 µl of PI in 5 ml of solution). The brain samples were subsequently incubated in labeling solution at 37°C with gentle shaking for 7–8 d. Next, the samples were incubated in washing solution containing 1.5% goat serum, 0.5% Triton X-100, and 0.05% sodium azide in 0.01 m PBS twice and overnight, and were additionally washed with 0.01 m PBS for 2 h three times at room temperature (RT). After completing the immunolabeling steps, tissue clearing was performed with the 3DISCO protocol (Ertürk et al., 2012a) by incubating the samples at RT with gentle shaking in gradient of tetrahydrofuran (THF; catalog #CP82.1, Carl Roth) water solutions (50 vol% THF, 70 vol% THF, 90 vol% THF, and 100 vol% THF overnight; and 100 vol% THF, 2 h for each step). Subsequently, the samples were immersed in dichloromethane (catalog #KK47.1, Carl Roth) for 1 h and finally in BABB (mixture of benzyl alcohol and benzyl benzoate in a 1:2 ratio; catalog #24122 and #W213802, Sigma-Aldrich) for at least 6 h until reaching full transparency. Aluminum foil was applied at all steps to avoid potential fluorescence quenching.
Whole brains of S24 and U87MG cohorts were optically cleared using the FluoClearBABB protocol (Breckwoldt et al., 2016), which is based on benzyl alcohol/benzyl benzoate clearing with a basic pH. Brains were dehydrated using analytical grade alcohol (t-butanol; Sigma-Aldrich) ranging from 30% to 100%. The clearing solution BABB consisted of benzyl alcohol (analytical grade; Merck) and benzyl benzoate (“purissimum P.A.” grade; Sigma-Aldrich) mixed in a 1:2 volume ratio. The pH levels of dehydration and clearing solutions were measured with an InLab Science electrode suited for organic solvents (METTLER TOLEDO) and adjusted with triethylamine (Sigma-Aldrich). Samples stored in BABB at 4°C in the dark were stable over time and showed no apparent decrease in fluorescence signal, allowing for reimaging.
Acquisition of light sheet microscopy datasets
Cleared brains of the GL261 cohort were scanned with a light sheet microscope (Luxendo LCS SPIM, Bruker) using a 4.0× objective lens (4×/340, 0.28 numerical aperture; XLFLUOR, Olympus) and combined lasers (excitation wavelength, 488, 561, and 642 nm; with respective filters and input power of 40 mW each). For image acquisition, an effective magnification of 4.4× with pixel size of 1.46 µm was used and tiling z-stack scans with a 6 µm step size were performed to cover the entire transparent brain samples. Total acquisition time was ∼60–70 min/brain with exposure times of 100 ms/slice. Images were stitched in Image Processor, exported as a series of TIF files, and used for further analysis. Brains of the S24 and the U87MG cohorts were scanned with a light sheet microscope (LaVision BioTec) using a 1.0× objective and a white light laser (SuperK EXTREME 80 mHz visible light, NKT Photonics). The wavelength spectrum ranged from 400 to 2400 nm (pixel size, 3.25 µm). The z-stacks with 5 µm step size and a total range of 1.5–2 cm were acquired for transversal coverage of the entire brain. Total acquisition time was ∼10 min prebrain with exposure times of 300 ms/slice. Images were exported as TIF files and used for further analysis.
Correlated MRI and LSM data analysis
Data processing and analysis was performed with open source 3D Slicer software [3D Slicer for Windows, versions 4.11.0 and 4.13.0 (www.slicer.org); Fedorov et al., 2012)]. Fractional anisotropy (FA) maps were calculated from DTI raw data using ParaVision 6 software (Bruker). Maps of T2 and T2* relaxation were calculated using an in-house MATLAB script (MATLAB R2020a for Windows, MathWorks). This included noise filtering with removal of all voxels exceeding 4 SDs of noise level of each dataset before fitting the data. MRI and LSM data were brought from the DICOM and TIF formats to the NIfTI format for further processing. To facilitate and speed up the processing pipeline, LSM data were downsampled with the bin shrink image filter of the Slicer Simple Filters module, applying a shrink factor of 5 × 5 × 5, which leads to the averaging of pixels within the 5 × 5 × 5 region into one pixel. The final resolutions of downsampled LSM data were 7.3 × 7.3 × 30 µm for the GL261 cohort and 3.25 × 3.25 × 25 µm for the U87MG and S24 cohorts, respectively. The Slicer Elastix extension was used to automatically register FA, T2, and T2* maps to the 3D T2w images of each animal. The resulting transform file was saved. A manual anatomic landmark-based registration was performed with the Fiducial registration wizard module of the Slicer IGT extension. Fiducials were placed by a neuroradiologist with >10 years of experience in small-animal MRI using anatomic structures that were unequivocally identifiable on each dataset (e.g., ventricles, blood vessels, or white matter tracts). A warping transform was created, and the resulting transform file was saved. Registration of T2w data to a mouse brain MRI atlas (Dorr et al., 2008) was performed in an analogous manner. Robustness of registration was tested by comparing the displacement vector field (DVF) of the final transform to DVFs of transforms established with subsets of the final fiducials in one animal. These subsets comprised every 10th, 5th, and 2nd of the final fiducials as well as only ventricle-, blood vessel-, and white matter-associated fiducials. Additionally, DVFs of output transforms generated with the Slicer Elastix extension with and without a manual initialization transform were compared with the DVF of the final landmark-based transform. All DVFs were generated in the Slicer transform module and compared using a dedicated in-house Python script computing the root mean squared error (rmse) for all the DVF grid points with respect to the optimal DVF. For each subset transform, the mean squared displacement with regard to the final landmark-based transform was calculated. Tumor segmentation was performed on T2w and LSM images in a semiautomated way using level tracing and threshold tools of the Slicer segment editor and was manually adjusted, if necessary. In addition, tumor subregions with densely packed tumor cells, tumor vasculature, and increased myeloid cell density were segmented on LSM data of the U87MG and the GL261 cohorts, respectively. The volume and extent of the segmentations were compared and mean FA values and T2 and T2* relaxation times were extracted, respectively. Moreover, mean FA values and T2 and T2* relaxation times were calculated for all anatomic regions provided in the mouse atlas. To minimize bias from differing spatial resolution of the atlas and the FA, T2, and T2* relaxation maps, only regions with a volume of at least 5 mm³ were included. This led to the inclusion of 20 of the 62 available anatomic atlas regions in the final analysis. Figure 1 provides an overview of the data processing and analysis pipeline.
Overview of the workflow. GL261, U87MG, and S24 glioma cells were orthotopically implanted. Multiparametric MRI data including T1-weighted and T2-weighted sequences, DTI, and T2 and T2* relaxometry were acquired. After scanning, brains were harvested, optically cleared, and imaged with LSM focusing on tumor cells, blood vessels, and myeloid cells. MRI and LSM datasets were coregistered to each other and to a reference mouse brain atlas, which differentiates 62 anatomic structures. Tumors and tumor subregions were segmented on T2-weighted and LSM images, respectively. Mean T2 and T2* relaxation times as well as FA were quantified, in the glioma, in histologically determined tumor subregions and in atlas-defined brain regions. To avoid partial volume effects, only atlas structures with a volume of at least 5 mm3 were included. This led to the inclusion of 20 of the 62 available anatomic atlas regions in the final analysis. Parts of the figure were created with Biorender.com.
Quantification of LSM data
Cx3cr1-EGFP-positive myeloid cells were quantified in the LSM data of GL261-inoculated mice. First, the segmentation of myeloid cells was performed using the Pixel Classification Workflow (Ilastik; Berg et al., 2019). The algorithm is based on a random forest classifier analysis. Cells and background were segmented by visual training until a visually satisfying result was obtained. Cells were segmented to obtain a single volume not connected to other volumes. The workflow yielded a segmentation map for each LSM dataset that was used as input for the subsequent Object Classification Workflow (Ilastik). Each connected volume was extracted. The workflow is based on a connected component analysis. The volumes were then filtered for size using a threshold of 3 voxels for nontumor regions and 4 voxels for tumor regions to reduce the false segmentation of noise or background. Because of the higher number of incorrectly segmented cells in the tumor region, a higher threshold was chosen. The information of segmented cells was then plotted in a map of the same size as the one used for manual registration using their coordinates to allow for multimodal comparison.
Based on visual inspection, anatomic areas were excluded from further analysis, which showed insufficient quality because of poor light sheet penetration or stray artifacts (consensus reading by two experienced readers, authors K.S. and M.O.B.). The frontal, parietotemporal and entorhinal cortex as well as the striatum, thalamus, hippocampus, and corpus callosum were selected for further analysis. To assess intratumoral heterogeneity, the entire LSM-defined tumor region was subdivided into tumor core and tumor periphery. For this, the “margin” function of the Slicer segment editor was used. “Tumor core” was generated by eroding the LSM-defined tumor label by 0–0.15 mm (corresponding to 0 × 0 × 0–20 × 20 × 5 voxels). For the tumor periphery, the LSM-defined tumor label was first dilated by 0–0.1 mm (corresponding to 0 × 0 × 0–14 × 14 × 3 voxels). Then, the tumor core label was subtracted from this dilated label to yield the final “tumor periphery” label. The extent of erosion and dilation was adapted for each animal based on visual inspection. To validate the segmentation data, volumes of 100 × 100 × 50 voxels were assessed by one reader (author K.S. with >10 years of experience in optical imaging) to visually determine the numbers of false-negative and false-positive segmented cells. Based on this, correction factors for the automated quantification were calculated. Further analysis of false-negative and false-positive rates as well as true cell count was performed independently for the tumor region and the nontumor region. Myeloid cell density of Cx3cr1-EGFP-positive cells/mm³ was calculated for all glioma and brain subregions.
Statistical analysis
Statistical analysis was performed with GraphPad Prism (version 9.3.1 for Windows, GraphPad Software). Median values of FA, T2, and T2* relaxation times of the LSM-defined tumors were compared with all anatomic atlas regions with Friedman's test followed by Dunn's multiple-comparisons test. Additionally, the median values of the quantitative MRI parameters in atlas-defined anatomic regions of the tumor-bearing and the contralateral hemispheres were compared with multiple paired Wilcoxon matched-pairs signed-rank test with correction for multiple comparisons using the Holm–Šídák method. Results from the automated quantification of myeloid cells were performed in an analogous manner. The comparison between tumor volumes derived from T2w and LSM segmentation in the three glioma models was performed using paired Wilcoxon matched-pairs signed-rank test with correction for multiple comparisons using the Holm–Šídák method as well. Comparisons between MRI parameters in T2w-segmented tumors and in LSM-defined subregions were realized with a Friedman's test followed by a Dunn's multiple-comparisons test. Differences between T2 and T2* relaxation times of GL261 and U87MG gliomas as well as FA between GL261 and S24 glioma were evaluated with unpaired Mann–Whitney tests. To compare the density of myeloid cells to MRI parameters in atlas-defined regions, Spearman nonparametric correlation was performed. Significance level was set to α = 0.05 and was corrected for multiple comparisons when appropriate. Multiplicity-adjusted p-values are reported wherever applicable.
Results
Landmark-based coregistration of MRI and LSM datasets is feasible
Coregistration of MRI and LSM datasets could not be accomplished using established automatic registration procedures because of the different contrasts and sizes of the images. Thus, we established an anatomic landmark-based approach with calculation of a warping transform based on manually placed fiducials on both datasets. This approach rendered optically satisfying results with good overlap of MRI and LSM data (Fig. 2A). To further test the robustness of this approach, we compared the mean squared displacement of transforms that were created with subsets of fiducials and automatic image registration to the final “ground truth” transform using fiducials in all identifiable structures (Fig. 2B). This comparison revealed that automatic coregistration rendered the worst result both qualitatively and in terms of the rmse of DVFs (Fig. 2B,C). Although the rmse values of the DVFs generated with subsets of fiducials placed on blood vessels, white matter tract-associated or ventricle-associated structures were in the range of only 0.28 mm (white matter tract-associated fiducials) to 0.19 mm (blood vessel-associated fiducials), and qualitative comparison showed suboptimal accuracy. For example, the use of only white matter tract-associated fiducials led to discrepancies along the cortical and midline structures of the brain. Application of only blood vessel-associated fiducials performed better along the cortical surface but did not align ventricles well (Fig. 2C). Thus, the best coregistration result was achieved by including fiducials placed on different anatomic structures.
Robustness of the coregistration. An anatomic landmark-based approach was used to coregister 3D T2w and corresponding LSM data of each animal. Fiducials were placed on anatomic structures, such as blood vessels, white matter tracts, and ventricles, that were unequivocally identifiable on both datasets. A, This rendered optically satisfying results with good overlap of MRI (in gray scale) and LSM data (in green scale). To evaluate this approach, coregistration was tested with transforms generated with subsets of the initial fiducial set. B, Overlays of T2w (in gray) and LSM data (in green) that were transformed with such subsets of fiducials. B, The transform generated with the Slicer Elastix module (left with purple frame) rendered the lowest accuracy of coregistration. The use of only white matter tract-associated fiducials led to discrepancies along the cortical and midline regions of the brain (middle with orange frame; the arrows mark areas lacking overlap). B, Application of only blood vessel-associated fiducials performed better along the cortical surface but did not align ventricles well (right with red frame; the arrows indicate suboptimal alignment of the ventricles). C, To further test the robustness of this approach, the rmse of the DVFs of transforms that were created with subsets of fiducials and the entire set of fiducials as ground truth were compared. The subsets of fiducials used are color coded, and the ground truth transform based on fiducials placed on all structures is marked with a star. C, Automatic coregistration using the Slicer Elastix module, without and with a manual initialization transform that approximately aligned both datasets, rendered the worst results not only qualitatively, but also in terms of the rmse (purple and light blue). B, C, Although the rmse of DVFs generated with subsets of fiducials placed on blood vessels (C, red), and white matter tract-associated (C, orange) or ventricle-associated structures (C, dark blue) were rather small, qualitative comparison showed suboptimal overlap (B).
Gray and white matter structures exhibit differences in MRI metrics, and sham injection has only local effects
To investigate the structural differences between gray and white matter structures, we compared FA, T2, and T2* relaxation times of atlas-defined anatomic regions in healthy control mice. T2 relaxation times of white and gray matter were comparable. Median T2* relaxation time was shorter in white matter such as the corpus callosum compared with gray matter like the cortex (median T2* relaxation time for corpus callosum, 29.34 ms; vs median T2* relaxation time for frontal cortex, 34.30 ms; Fig. 3A). T2* relaxation maps of the entorhinal cortex were prone to susceptibility artifacts caused by air-filled structures of the temporal bone. Therefore, the entorhinal cortex was excluded from further analyses of regional T2* relaxation time. FA of white matter structures was overall higher than in gray matter structures (median FA corpus callosum, 0.323; vs median FA frontal cortex, 0.188; Fig. 3A). There were no significant differences between the hemispheres in any of the parameters investigated. Sham injection of PBS in the right striatum led to a visible hypointense tract surrounded by a small hyperintense rim on T2w images, which can be ascribed to blood products in and edema around the injection site. This was reflected in a significantly lower T2* relaxation time of the needle tract compared with most brain regions (Fig. 3B). Median T2 relaxation time of the needle tract was also shorter than in some regions, while median FA of the tract was comparable to normal brain parenchyma (Fig. 3B). However, the presence of the needle tract did not affect global MRI metrics of the brain region it was located in, as median FA, T2, and T2* relaxation times of the right PBS-injected and the left healthy hemisphere remained similar.
Atlas-based analysis of healthy brain and of the effects of PBS sham injection. A, Comparison of T2 and T2* relaxation times as well as FA between brain regions in the left and right hemisphere of healthy C57BL/6 mice. As expected, none of the parameters differ significantly between the hemispheres. FA of white matter tracts such as the corpus callosum is higher than in gray matter structures. B, We investigated the effects of PBS sham injection on regional MRI parameters. The injection site could be identified as a hypointense tract on T2w images and was surrounded by a narrow T2w hyperintense rim. B, This qualitative impression is reflected by significantly shorter T2* and T2 relaxation times of the PBS injection tract compared with brain regions. However, the presence of the injection site had no impact on global regional MRI features, as the MRI parameters of the right PBS-injected hemisphere and the left healthy hemisphere remained comparable. PBS injection had no effect on local or regional FA. The asterisks indicate the level of significance of multiplicity adjusted p-values (α = 0.5).
GL261 glioma causes regional differences in quantitative MRI metrics, its subregions exhibit different imaging characteristics and T2 relaxation time correlates with myeloid cell density
To apply our pipeline in a disease context, we conducted a comprehensive multiparametric MRI protocol in the syngeneic GL261 glioma model (Fig. 4A) as a commonly used immunocompetent murine model of glioblastoma (Tritz et al., 2021). When comparing the quantitative MRI parameters of the tumors segmented on T2w and LSM images, respectively, to the anatomic regions, FA and T2 relaxation time were the parameters that best distinguished the tumor from atlas-defined regions. Median FA of the tumor was lower than median FA of the corpus callosum, hypothalamus, midbrain, and pons. The GL261 glioma had a longer T2 relaxation time than the contralateral striatum, the corpus callosum, and the thalamus, while T2* relaxation time was shorter than in the contralateral frontal lobe (Fig. 4B–D). Further analysis of quantitative MRI parameters did not reveal significant differences when comparing the atlas-defined anatomic regions of the hemisphere, in which the tumors were implanted (“inoculated” in the following) with the corresponding regions in the contralateral hemisphere. However, there was a trend toward longer T2 relaxation and higher FA in tumor-bearing regions (Fig. 4B–D, Table 2).
Comparison of anatomic brain regions with suspected tumor on T2w MR images and confirmed tumor on LSM data for the three glioma models
Atlas-based analysis of the GL261 cohort. A, Exemplary images from an animal with a GL261 glioma. Of note, MR images are not displayed in radiologic convention to match the LSM data (i.e., the right hemisphere is on the right side of the image). T2, T2*, and FA maps with an overlay of atlas labels displaying the anatomic regions analyzed (outlined in different colors) and the segmented tumor (pink area) are presented from top to bottom. The striatum (light green) and the frontal cortex (light blue) contained the largest tumor burden. B–D, T2 relaxation time (B) and T2* relaxation time (C) as well as FA (D) were quantified in atlas-defined anatomic regions and in GL261 tumors segmented on LSM datasets (red bar and dots). The hemisphere in which the tumor was implanted is displayed in gray (bars and dots) and the contralateral hemisphere is shown in white (bars and dots). Unpaired brain regions are displayed in dark gray (bars and dots). Regions with LSM-confirmed presence of tumor cells are marked by red stripes. Median values in the LSM-confirmed tumor region were compared with all atlas-defined structures. Then, median values of the quantitative MRI parameters of the inoculated hemisphere were compared with the contralateral side. Unpaired brain regions (medulla, midbrain, and pons) were not included in this analysis. B, T2 relaxation time of GL261 glioma was longer than in several brain regions. C, D, Moreover, T2* relaxation time (C) and FA (D) of the tumor differed from anatomically defined regions. T2 relaxation time and FA distinguished best between glioma and brain tissue. The asterisks indicate the level of significance of multiplicity adjusted p-values. The α value was set to 0.5.
LSM analysis was used to identify subregions in a GL261 glioma. GL261 cells were not evenly distributed within the bulk tumor. Areas with densely packed tumor cells were visually identified and segmented based on higher fluorescence compared with the surrounding tumor with interspersed, lower tumor cell density (Fig. 5A). These highly cellular regions and areas with accumulation of Cx3cr1-positive myeloid cells were distinguished from the remaining tumor. Visual inspection of coregistered T2w and LSM data showed that the areas of T2w-suspected and LSM-confirmed tumor highly overlapped. However, the myeloid cell accumulation extended the main tumor region and was also present in the adjacent brain parenchyma (Fig. 5A). Median FA, T2, and T2* relaxation times were analyzed in these subregions and compared with median values derived from the entire tumor as segmented on T2w images. Median FA of subregions with dense tumor cells and accumulated myeloid cells were higher than in the entire tumor (Fig. 5B). The median T2 relaxation time of areas with high myeloid cell density was significantly shorter than in the T2w segmented entire tumor (Fig. 5B).
Comparison of GL261 subregions. A, Coregistered MRI and LSM data of an animal with a GL261 glioma. A, The tumor was segmented on T2w images (A, top left, light red). Scale bar, 1 mm. The inlet shows a magnification of the tumor, which grows bulky and is mainly T2w hyperintense compared with the surrounding tissue. Scale bar, 1 mm. LSM datasets were used to define subregions of a GL261 glioma: areas with densely packed tumor cells (A, top middle; green; scale bar, 500 µm) were differentiated from the remaining tumor (A, top middle, yellow; scale bar, 500 µm) and from areas with accumulated Cx3cr1-EGFP-positive myeloid cells (microglia, monocytes, dendritic cells, and NK cells; A, top right, pink; scale bar, 500 µm). A, Magnification of the T2w image (left inlet) reveals a rather homogeneously T2w-hyperintense tumor center with a T2w-hypointense rim. Scale bar, 1 mm. A, Magnification of the LSM images reveals that tumor cells are not evenly distributed, but rather exhibit a radial formation with interspersed areas of lesser cell density (middle inlet). Scale bar, 100 µm. A, Myeloid cells accumulate along the tumor borders, but can also be found in the adjacent tissue (right inlet). Scale bar, 100 µm. Overlays of these regions on maps of T2 relaxation time and FA are shown below. The areas of T2w-suspected (light red) and LSM-confirmed tumor (yellow/green) overlap, while myeloid cell accumulation (pink) extends beyond the main tumor bulk. B, When comparing quantitative MRI parameters between these regions, significant differences were observable between mean values of the T2w-segmented tumor and LSM-defined subregions. C, The density of Cx3cr1-EGFP-positive myeloid cells was automatically quantified. On the left, an exemplary segmentation output is shown (white, Cx3cr1-EGFP-positive myeloid cells; downsampled voxel size, 7.3 × 7.3 × 30 µm). Scale bar, 500 µm. Millions of myeloid cells, in this case 2,096,833 cells, were automatically segmented, which showcases the power of our segmentation approach. By this, cell quantities in the tumor and in anatomic brain regions were made accessible for further analyses. On the right, a heat map displaying the number of myeloid cells per voxel of a downsampled LSM dataset is shown (voxel size, 255.5 × 255.5 × 1050 µm). Scale bar, 500 µm. The highest number of myeloid cells per voxel was observable in the tumor periphery. D, Myeloid cell density per mm³ was calculated and compared between brain regions. The tumor periphery displayed the highest myeloid cell density, significantly exceeding most non-tumor-bearing regions. D, Myeloid cell density was lowest in the thalamus. E, We additionally correlated regional myeloid cell density to T2 and T2* relaxation time as well as FA and observed a significant positive correlation between myeloid cell density and T2 relaxation time. Asterisks indicate the level of significance of multiplicity-adjusted p-values.
Moreover, we automatically quantified the density of Cx3cr1-positive myeloid cells in the tumor core, the tumor periphery, and in a subset of atlas-defined anatomic regions. Myeloid cell density was highest in the tumor periphery (Fig. 5C) and showed significant differences when compared with the contralateral striatum, corpus callosum, entorhinal cortex, as well as the thalamus in both hemispheres (Fig. 5D). Again, there was a trend of higher myeloid cell content in tumor-bearing regions (Fig. 5D). Interestingly, there was a significant positive correlation between T2 relaxation time and myeloid cell density (Spearman's r = 0.76, p = 0.0368; Fig. 5E), whereas FA and T2* relaxation time did not correlate with myeloid cell density.
U87MG gliomas have a different MRI signature than GL261
U87MG gliomas grow bulky and exhibit intratumoral heterogeneity and extensive neoangiogenesis (Breckwoldt et al., 2019). As the presence of pathologic vasculature is more likely represented in relaxation time rather than FA, the MRI protocol was focused on T2 and T2* relaxometry. In line with the GL261 glioma, the mean T2 relaxation time was also the best parameter to distinguish tumor and all other brain regions in a U87MG glioma (Fig. 6A). There were no significant differences between T2* relaxation time in U87MG gliomas and atlas-defined subregions (Fig. 6B). When comparing atlas-defined anatomic regions of the inoculated to the contralateral hemisphere, there was a trend toward lower T2* relaxation time in tumor-bearing regions (Fig. 6A,B, Table 2). In contrast to the GL261 model, imaging characteristics in subregions of a U87MG glioma identified on LSM data were similar to those in the whole tumor (Fig. 6A,B).
Atlas-based and subregional analyses of the U87MG and S24 cohorts. A–C, Analysis of T2 relaxation time (A, left) and T2* relaxation time (B, left) as well as FA (C) in atlas-defined anatomic regions and in U87MG tumors (A, B) and S24 tumors (C) segmented on LSM datasets (red bars and dots). The hemisphere in which the tumor was implanted is displayed in gray (bars and dots), and the contralateral hemisphere is shown in white (bars and dots), while unpaired regions are shown in dark gray (bars and dots). Regions with an LSM-confirmed presence of tumor cells are marked by red stripes. Quantitative MRI parameters of the LSM-segmented tumor were compared with all atlas regions. Additionally, regions of the inoculated hemisphere were compared with the contralateral side. A, B, Additionally, mean T2 relaxation time (A, right) and T2* relaxation time (B, right) of subregions defined on LSM datasets (green, densely packed tumor cells; red, tumor vessels; yellow, remaining tumor) were compared with each other and to the entire tumor, but T2 and T2* relaxation times were similar in these regions. The asterisks indicate the level of significance of multiplicity adjusted p-values (α = 0.5).
In contrast to GL261 and U87MG, which show encapsulated growth, the third used the S24 glioma model, which forms large networks of connected tumor cells and diffusely infiltrates the brain parenchyma using white matter tracts as guiding structures (Osswald et al., 2015; Venkataramani et al., 2019). As tumor cell invasion most likely affects tissue microstructure and should therefore be detectable in diffusion parameters, the MRI protocol included DTI, and the analysis was focused on FA. In contrast to the other models, S24 glioma cells diffusely infiltrate the brain rather than grow as tumor bulks. Neoangiogenesis was not detected based on lectin staining of the microvasculature (data not shown). The comparison of atlas-defined brain regions showed no significant differences between hemispheres or between the tumor and anatomic regions (Fig. 6C). The striatum of the inoculated hemisphere was the brain region that carried the largest tumor burden and tended to have a lower FA than the contralateral counterpart (Fig. 6C, Table 2).
Overall, the U87MG glioma had a different MRI signature than the GL261 glioma. Both T2 and T2* relaxation times of U87MG tumors were significantly lower (median T2 relaxation time: 48.37 vs 51.08 ms; p = 0.0317; median T2* relaxation time: 14.48 vs 28.76 ms; p = 0.0079; Fig. 7A). Mean FA values of GL261 and S24 gliomas were not significantly different.
Comparison of MRI-characteristics and tumor volume of the three glioma models. A, GL261 gliomas (white bars and dots) and U87MG gliomas (light gray bars and squares) present with significantly different T2 and T2* relaxation times (left and middle). FA values of GL261 and S24 gliomas (dark gray bar and triangles) were similar. B, T2w images with segmented tumors and coregistered slices from LSM data displaying the confirmed tumor region are depicted for each glioma type. Scale bars, 500 µm. Additionally, overlays of coregistered T2w and LSM images (in red for better visibility) are shown. Tumor volumes derived from T2w-based (pink bars and dots) and LSM-based (red bars and dots) segmentation were compared. B, C, The volumes of T2w-segmented and LSM-segmented tumors overlap in bulky growing gliomas such as GL261 and U87MG (B, top and middle row, C). In these tumor types, T2w suspect areas can be used as a surrogate for tumor presence. In contrast, the T2w-suspicious and LSM-confirmed tumor areas do not completely overlap in invasively growing S24 glioma, leading to a significant difference between the volumes derived from T2w and LSM segmentation (B, bottom row, C). Thus, T2w segmentation of S24 glioma does not necessarily only reflect tumor-bearing parenchyma, but may also contain edema, for example.
The extent of S24 glioma infiltration cannot be reliably detected by conventional MRI
Delineation of the S24 glioma on T2w images was more difficult than for the GL261 and U87MG gliomas, which had an easily identifiable, encapsulated, and bulky growth pattern (Fig. 7B). In contrast, the S24 glioma presented with subtle T2w hyperintensity and mass effect (Fig. 7B), but borders remained ill defined in all acquired MRI sequences. Therefore, T2w segmentation of the S24 glioma included the entire area of faint T2w hyperintensity, which was suspected to be tumor infiltrated based on visual MRI inspection. In contrast, LSM-based segmentation of regions containing fluorescent S24 cells was unambiguous and was used as a cellular ground truth of tumor extent. Interestingly, the comparison of tumor volumes derived from T2w-based and LSM-based segmentation showed significant differences between T2w-suspected and LSM-confirmed tumor regions (T2w-suspected volume, 105.9 mm³; vs LSM-confirmed volume, 36.3 mm³; p = 0.0081; Fig. 7C), whereas T2w-suspected and LSM-confirmed tumor volumes of GL261 and U87MG gliomas were almost identical (Fig. 7C). Thus, we identify different MRI and LSM signatures in these brain tumor models.
Discussion
MRI is paramount for the initial diagnosis and treatment monitoring of a variety of neurologic diseases, including brain tumors. However, MRI has shortcomings regarding sensitivity, specificity, and resolution. Advancement in MRI technology (field strength, coils, and sequence development) has improved in particular sensitivity and in parts also specificity. To further understand the molecular and cellular basis of MR signals, an important limitation remains the lack of ground truth information for validating multiparametric MRI, thus correlating key cellular and molecular features of the disease with in vivo imaging findings (Hempel et al., 2019). To overcome these limitations, correlative histologic analyses of glioma have been performed in depth (Zaccagna et al., 2019; Bobholz et al., 2020). However, similar to various other studies, for example of the developing brain (Breu et al., 2019), these studies are often restricted to 2D analysis of single tissue sections or regions of interests, thus not reflecting the highly complex 3D anatomy of the brain.
Tissue clearing and light microscopy have developed as an entire field to drive discovery in neuroscience (Ueda et al., 2020). LSM offers 3D information of the entire brain with submicrometer resolution, and basic neuroscience questions have been investigated in depth illustrating, for example, neuronal projections for dissecting specified neuronal circuits in whole organs and even organisms. LSM has also been applied to models of pathology and neuropathology (Ertürk et al., 2012b), has been used in cancer research (Almagro et al., 2021), and lately has entered human pathology (Zhao et al., 2020). Previous work also used information gained from LSM to analyze and interpret MRI data. For example, a multimodality and multiscale imaging and visualization pipeline for vascular systems was developed by Bhargava et al. (2022). While a similar fiducial-based landmark registration of images was used, MRI was performed ex vivo and only 1.5 mm sections of the mouse brain were optically cleared and coregistered, contrasting this work to our approach. Another work used LSM data as substrates for computational modeling, which was guided by in vivo MRI data (d'Esposito et al., 2018). While information from both LSM and MRI was used to predict drug uptake and treatment responses in murine tumor models, the authors did not coregister or correlate MRI and LSM directly.
We have previously shown how MRI and LSM (ultramicroscopy) datasets can be acquired side by side to inform MRI on cellular correlates as shown by LSM (Breckwoldt et al., 2016, 2019). In the current work, we further improve this pipeline by voxel-to-voxel coregistration of both datasets and further coregistration with an anatomic atlas for standardized quantitative imaging analysis of MRI and LSM datasets. Moreover, we established an automated segmentation and quantification pipeline for the entire myeloid compartment in GL261 glioma-bearing animals. This allowed single-cell segmentation of Cx3cr1-positive myeloid cells, thus enabling the identification and quantification of millions of single cells. Our method made the complex three-dimensional LSM data with an average size of ∼100 GB and encompassing a median number of 1.9 × 106 macrophages/microglia in the analyzed regions accessible for direct correlation with MRI metrics.
Atlas-based analyses revealed differences in quantitative MRI and LSM parameters between tumors and gray and white matter subregions in the inoculated and contralateral hemisphere in all three glioma models. Moreover, MRI parameters and the myeloid cell density differed in structures that were infiltrated with tumor cells as confirmed by LSM compared with their contralateral counterparts. This can be interpreted as the effects of glioma and their microenvironment on global and regional imaging features. The lack of statistical significance when comparing subregions of both hemispheres is likely because of the fact that the entire anatomic region was not occupied by a glioma. This rigorous anatomic and unbiased analysis of quantitative MRI metrics can be particularly useful in tumor models with diffuse infiltration of the adjacent brain parenchyma like S24, as these are difficult to delineate on conventional MRI sequences and thus hamper tumor segmentation on MRI. The combination of atlas-based analyses and coregistered LSM images allows the unambiguous identification of tumor-bearing regions to facilitate the interpretation of MRI results as brain regions have inherent differences in their MRI signal characteristics (white vs gray matter, fiber tracts vs deep nuclei; Qin et al., 2013).
The coregistration of MRI and LSM datasets enables the analysis of cellularly resolved, distinct subregions. Our results show that such subregions can differ in their FA and T2 relaxation times from the entire tumor. Moreover, we identified a correlation of regional myeloid cell density and T2 relaxation time. Therefore, subregional analyses can contribute to an in-depth characterization of a glioma. Our approach is also compatible with advanced radiomic-based and artificial intelligence-based analyses and future work should aim to identify distinct MRI signatures of such subregions and to translate these findings to clinical imaging of a human glioma. Ultimately, this could aid treatment planning and monitoring of a glioma.
T2w-based segmentation of GL261 and U87MG gliomas was in good agreement with LSM-based segmentation. In contrast, tumor volumes derived from T2w and LSM segmentations were markedly different in the S24 model. The area in which tumor was suspected on T2w images was larger than the LSM-confirmed tumor region. A possible explanation is that the T2w segmentation included not only tumor mass but also associated edema. The differentiation of tumor borders is a well appreciated clinical problem in glioma imaging and treatment. “Peritumoral” regions of a human glioblastoma consist of infiltrating tumor cells, astrogliosis, vasogenic edema, and an active tumor immune microenvironment (Akbari et al., 2016), which cannot be faithfully differentiated based on conventional T2w or fluid-attenuated inversion recovery sequences. The same phenomenon is conceivable for the S24 glioma model. Another possibility is that some fluorescent signal of S24:td-tomato was lost during the clearing process, so that the LSM-confirmed tumor region was underestimated. A combination of both issues is also possible, and we appreciate that this discrepancy warrants further investigation. Meanwhile, T2w hyperintensity as a surrogate for infiltrative glioma should be used with caution.
Although we tested our approach in three different glioma models, further research is needed to prove generalizability. Also, the use of endogenous fluorescent proteins such as td-tomato or GFP in combination with exogenous fluorescent dye conjugates can lead to varying preservation of fluorescent signals during the clearing process. Moreover, the direct translation of our work to human clinical studies is difficult. First, our work is based on correlative MRI and light sheet microscopy data, which are not readily available in the clinical setting and require a surgically resected or autopsy specimen for LSM. We have previously shown that correlative clearing and high-field MRI of human samples is possible (Breckwoldt et al., 2019). These tissue specimens are, however, currently limited to 1–2 cm³. With the recent development of whole-organ clearing (Zhao et al., 2020), our pipeline could in principle also be performed in human organs postautopsy. We speculate that our fiducial-based coregistration pipeline should also work for such datasets acquired at 3 T, but this needs to be investigated in future studies.
In summary, the direct correlation of MRI and LSM data allows for an in-depth characterization of glioma models and can be used to decipher the structural, cellular, and molecular foundation of tumoral MRI signals. The proposed analysis pipeline can be applied to other preclinical brain tumor or neurologic disease models. We envision that imaging signatures derived from such analyses could also ultimately improve the diagnosis and management of patients with brain tumors as they enhance the interpretation of MRI.
Footnotes
K.S. is supported by the Olympia-Morata program of the Medical Faculty, University of Heidelberg, and acknowledges funding from the Daimler and Benz Foundation (Grant 32-97/21). M.O.B. acknowledges funding from the Emmy Noether Program of the German Research Foundation (DFG; Grant BR 6153/1-1) and the Else Kröner-Fresenius Stiftung (Grants 2017-A25 and 2019_EKMS.23). W.W., T.B., M.P., F.W., M.B., and M.O.B. were funded by the UNITE Glioblastoma (Project-ID 404521405; Grant SFB 1389; work packages A03, B03, C02, and C03). C.P. is funded by the DFG (Grant 451894423). V.V. received financial support from the DFG (Grant VE1373/2-1), the Else Kröner-Fresenius Stiftung (Grant 2020-EKEA.135), and the University of Heidelberg (Physician Scientist-Program and Krebs und Scharlachstiftung). We thank Dr. Gergely Solecki (DKFZ Heidelberg) for help with tumor cell xenografting.
The authors declare no competing financial interests.
- Correspondence should be addressed to Michael O. Breckwoldt at michael.breckwoldt{at}med.uni-heidelberg.de or Katharina Schregel at katharina.schregel{at}med.uni-heidelberg.de
