Multidimensional MRI for Characterization of Subtle Axonal Injury Accelerated Using an Adaptive Nonlocal Multispectral Filter

1 Introduction

Water molecules within biological tissues interact with their local chemical environment via nuclear relaxation processes and follow diffusion patterns trajectories that are governed by the local tissue density and geometry. Using a combination of magnetic field profiles to probe these mechanisms, magnetic resonance (MR) provides exquisite sensitivity to both the chemical composition, through relaxation parameters, and microstructure, through diffusion parameters, of biological tissues.

One fundamental obstacle for using MRI to characterize tissue heterogeneity is the averaging that occurs across the image volume elements, known as voxels (i.e., pixels with thickness). Voxel-averaged images can only provide macroscopic information with respect to the voxel size, which is typically ∼1–3 mm³. In a mammalian brain, an individual voxel contains multiple chemical and physical microenvironments such as axons, neurons, glia, myelin, and cerebrospinal fluid. Many biological processes-of-interest take place at a microscopic scale that only affects a small portion of any given voxel, which therefore makes them undetectable using conventional voxel-averaged MRI methods. The inability to separate normal and pathological tissue within a voxel is a major contributor to the insensitivity and ensuing non-specificity of conventional MRI methods in detecting abnormal cellular processes.

By simultaneously encoding multiple MR “dimensions”, such as relaxation times (T₁ and T₂) [1] and diffusion [2, 3], multidimensional distributions of those MR parameters can provide fingerprints of various chemical and physical microenvironments within the volume-of-interest, which can be traced back to specific materials and cellular components. If combined with imaging [4], multidimensional MRI has the potential to overcome the voxel-averaging limitation by accomplishes two fundamental goals: 1) it provides unique intra-voxel distributions instead of an average over the whole voxel; this allows identification of multiple components within a given voxel [5–7], while 2) the multiplicity of dimensions inherently facilitates their disentanglement; this allows higher accuracy and precision in derived quantitative values [8–11].

Although traditionally multidimensional MR experiments required many repeated acquisitions and therefore have imposed serious time constraints [12], acquisition strategy [13, 14], computational [3, 6, 15, 16], and pulse design [17, 18] technological breakthroughs have significantly reduced the data burden and positioned multidimensional MRI as a powerful emerging imaging modality for studying biological media. Despite of these advances, wide-spread clinical translation still presents challenges, in particular, due to relatively low signal-to-noise ratio (SNR) and the ensuing increased data amount requirement. To address that, we report the use of a nonlocal estimation of multispectral magnitudes (NESMA) filter [19] on multidimensional MRI data to perform noise reduction for reliable parameter determination and further data reduction. To date, NESMA has been successfully used to improve determination of myelin water fraction from multi-spin-echo MR images [20], or cerebral blood flow from arterial spin labeling MR images [21].

We chose to focus and reprocess a subset of data from our recent study that showed multidimensional MRI can uncover subtle axonal injury patterns in the human brain, otherwise inaccessible using conventional quantitative MRI techniques such as diffusion tensor imaging (DTI), T₁ or T₂ maps [22]. The study investigated brain samples derived from human subjects who had sustained traumatic brain injury (TBI) and control brain donors using MRI, followed by co-registered histopathology that included amyloid precursor protein (APP) immunoreactivity to define axonal injury severity [23]. Abnormal multidimensional T₁-T₂, mean diffusivity-T₂ (MD-T₂), and MD-T₁ spectral signatures that were strongly correlated with injured voxels were identified and used to generate axonal injury biomarker maps [22]. Here we study the effect of applying a multispectral nonlocal filter on three representative cases (a control and two TBI cases), with the main goal of evaluating the performance of NESMA on reduced datasets to limit the total acquisition time required for reliable multidimensional MRI characterization of brain tissue. The co-registered APP histology images serve as a “ground truth” reference, thus providing a unique opportunity to quantitatively evaluate to what extent the accuracy of the injury biomarkers maps is preserved under substantial data reduction.

2 Methods

2.1 Donors Specimens Employed in the Present Study

We evaluated autopsy-derived brain specimens from two different human brain collections. Formalin-fixed portions of approximately 20 × 20 × 10 mm³ of the Corpus Callosum (CC) were obtained from one military subject from the DoD/USU Brain Tissue Repository and Neuropathology Program (https://www.researchbraininjury.org, Bethesda, MD; Subject 1), and two civilian subjects enrolled in the Transforming Research and Clinical Knowledge in Traumatic Brain Injury study (TRACK-TBI; https://tracktbi.ucsf.edu/transforming-research-and-clinical-knowledge-tbi) (Subjects 2 and 3). For each case, a next–of-kin or legal representative provided a written consent for donation of the brain for use in research. The brain tissues used have undergone procedures for donation of the tissue, its storage, and use of available clinical information that have been approved by the USU Institutional Review Board (IRB) prior to the initiation of the study. All experiments were performed in accordance with current federal, state, DoD, and NIH guidelines and regulations for postmortem analysis.

Subject 1 was a 44 years old male with no known TBI history and postmortem APP-negative histopathology. Subject 2 was a 60 year old male that died as a result of a intraparenchymal hemorrhage following a motor vehicle accident. Subject 3 was a 49 year old male that died as a result of intraparenchymal and subarachnoid hemorrhages following a fall.

2.2 MRI Acquisition

Prior to MRI scanning, each formalin-fixed brain specimen was transferred to a phosphate-buffered saline (PBS) filled container for 12 days to ensure that any residual fixative was removed from the tissue. The specimen was then placed in a 25 mm tube, and immersed in perfluoropolyether (Fomblin LC/8, Solvay Solexis, Italy), a proton free fluid void of a proton-MRI signal. Specimens were imaged using a 7 T Bruker vertical bore MRI scanner equipped with a microimaging probe and a 25 mm quadrupole RF coil.

Multidimensional data were acquired using a 3D echo planar imaging (EPI) sequence with a total of 56 and 302 images for T₁-T₂ and MD-T₂, respectively, and with 300 μm isotropic spatial resolution, which resulted in respective acquisition times of 4.5 and 17.8 h. To test the feasibility of data reduction using NESMA we derived reduced datasets by sub-sampling the full datasets. The total number of T₁-T₂ images was reduced from 56 to 36 (35.7% decrease), while the total number of MD-T₂ images was reduced from 302 to 122 (59.6% decrease). Further details can be found in the Supplementary Material.

The SNR was always maintained above 100 (defined as the ratio between the average unattenuated signal intensity within a tissue region of interest, and the standard deviation of the signal intensity within the background). The sample temperature was set at 16.8°C.

2.3 Multidimensional MRI Processing

Here we implemented a marginally-constrained, ℓ₂-regularized, nonnegative least square optimization to compute the multidimensional distribution in each voxel, as previously described [8, 24]. It is a well-tested approach that had been proved robust and reliable [2, 14, 25–29], which in this study had resulted in two types of distributions in each voxel: T₁-T₂ and MD-T₂. The 2D T₁-T₂ and MD-T₂ distributions were evaluated on 50 × 50 logarithmically sampled grids using a previously described algorithm [13]. The range for T₁ was 1–10,000 m, the range for T₂ was 1–500 m, and the range for MD was 0.0001–5 μ m²/ms.

If one considers the multidimensional distributions as spectra, it is possible to use them to generate maps of specific spectral components by means of integration over a predefined parameter range generally associated with a spectral peak. The integral value is a number between 0 and 1, representing a certain spectral component (SC) in a given multidimensional distribution, which can be computed in each voxel to generate an image of that specific SC [30]. Here we apply a recently proposed unsupervised algorithm to identify the injury-associated spectral information [22], and generate injury biomarker maps that closely follow APP histopathology.

2.4 The Nonlocal Estimation of Multispectral Magnitudes (NESMA) Filter

For each sample, the multidimensional distributions were derived from the original multidimensional data as well as from data denoised using the NESMA filter to improve accuracy and precision in derived distributions.

We consider K multidimensional images defined on a discrete grid describing the 3D spatial domain spanned by the image. The underlying idea of quantitative filters is to reduce noise by replacing the intensity of a given voxel by an unbiased estimate of its underlying amplitude. This requires selection of voxels that are likely to come from similar tissue. The NESMA filter restores the amplitude, A, of an index voxel, i, based on intensities of M preselected voxels with similar multispectral signal patterns through:

$A_k(i)=\frac{1}{M}\sum _j^M{S}_k(j),(1)$

where S_k(j) is the measured amplitude in voxel j of frame k. M is the total number of similar voxels defined using the relative Manhattan distance (RMD) between voxel intensities as

$RMD(i,j)=100\times \frac{{\sum }_{k=1}^K\left|S_k(i)-S_k(j)\right|}{{\sum }_{k=1}^K{S}_k(i)}.(2)$

The RMD was calculated between the index voxel i and all voxels belonging to a relatively large search window of size R, centered around the index voxel i, in which emission and reception B₁ fields and noise standard deviation (SD) were assumed to be approximately constant. The size of the window must be sufficiently large to ensure inclusion of an adequate number of similar voxels, and sufficiently restricted to ensure that the transmission and reception B₁ fields and noise SD are approximately constant within the window. In this work, we used a relatively conservative window size to avoid introducing bias in the estimated amplitudes. The size of the search window, R, was fixed at 11 × 11 × 11 voxels. Voxels with RMD <5% were considered similar to the index voxel.

2.5 Histopathology

After MRI scanning, each CC tissue block was transferred for histopathological processing. Tissue blocks from each brain specimen were processed using an automated tissue processor (ASP 6025, Leica Biosystems, Nussloch, Germany). After tissue processing, each tissue block was embedded in paraffin and cut in a series of 5 μm-thick consecutive sections on which immunohistochemistry for anti-amyloid precursor protein (APP) was performed (DS9800, Leica Biosystems, Buffalo Grove, IL). Further details can be found in [22].

2.6 Quantification of Axonal Damage

Images of APP stained sections were digitized using an Aperio whole slide scanning scanner system (Leica Biosystems, Richmond, IL) at ×20 magnification. The following steps, all implemented using MATLAB (The Mathworks, Natick, MA), were taken to allow for a quantitative analysis of the APP images. First, the images were transformed into a common, normalized space to enable improved quantitative analysis [31]. Then, the normalized images were deconvolved to unmix the primary (APP) and secondary (hematoxylin and eosin, H&E) stains, and background to three separate channels [32]. Once an APP-only image was obtained, a final thresholding step was taken to exclude non-specific staining and to allow for a subsequent % area calculation.

From each tissue section, based on APP staining, traumatic axonal injury (TAI) lesions were identified by an experienced neuropathologist (DI) as white matter (WM) areas with swollen axonal varicosities, axonal bulbs, or distorted axons. Accordingly, regions of interest (ROI) of normal-appearing WM and TAI lesions were manually defined. Additionally, gray matter (GM) ROIs from adjacent cingulate cortex were defined in all sections. Twelve ROIs, covering together an average 81 mm² of tissue, were identified per tissue section. After extracting the ROIs, APP density was expressed as the percentage of total area within the ROI in the binary deconvolved APP image. In total, 36 ROIs from three subjects were included in this study.

3 Results

3.1 Axonal Injury Spectral Signatures Are Preserved After Filtering

We first investigated the spatially-resolved subvoxel T₁-T₂ and MD-T₂ spectral components to assess the effect of NESMA on the derived voxelwise spectra. To do that, it is useful to summarize the 4D information, which consists of 2D images with 50 × 50 spectra in each voxel, as arrays of images with varying subvoxel T₁, T₂, and MD values. To make them more readable, the 50 × 50 spectra were sub-sampled on a 10 × 10 grid. These maps are shown in Figures 1–3 for all three Subjects. Corresponding histological APP images (co-registered with the MRI) are shown on the left panel of Figure 4, with red color indicating abnormal APP accumulation.

FIGURE 1

FIGURE 1. Maps of 2D probability density functions (i.e., 2D normalized spectra) from Subject 1 (control) of (A) unfiltered and (B) filtered subvoxel T₁-T₂ values reconstructed on a 10 × 10 grid of subvoxel T₁ values (horizontal axes) and subvoxel T₂ values (vertical axes), and maps of (C) unfiltered and (D) filtered subvoxel MD-T₂ values reconstructed on a 10 × 10 grid of subvoxel MD values (horizontal axes) and subvoxel T₂ values (vertical axes).

Starting with the control case (Subject 1), the spatially-resolved subvoxel T₁-T₂ and MD-T₂ spectral components are shown in Figure 1. The left column shows the results from the unfiltered data (T₁-T₂ and MD-T₂ in Figures 1A, C, respectively), while the right column shows the results from the filtered data (T₁-T₂ and MD-T₂ in Figures 1B, D, respectively). The maps revealed signal components that were spatially consistent with specific tissue types such as white matter and gray matter.

The spatially-resolved subvoxel T₁-T₂ and MD-T₂ spectral components from the first TBI case (Subject 2) are shown in Figure 2. As before, the left column shows the results from the unfiltered data (T₁-T₂ and MD-T₂ in Figures 2A,C, respectively), and the right column shows the results from the filtered data (T₁-T₂ and MD-T₂ in Figures 2B,D, respectively). Similarly to the control case, here too the maps revealed signal components that were spatially consistent with specific tissue types.

FIGURE 2

FIGURE 2. Maps of 2D probability density functions (i.e., 2D normalized spectra) from Subject 2 (TBI) of (A) unfiltered and (B) filtered subvoxel T₁-T₂ values reconstructed on a 10 × 10 grid of subvoxel T₁ values (horizontal axes) and subvoxel T₂ values (vertical axes), and maps of (C) unfiltered and (D) filtered subvoxel MD-T₂ values reconstructed on a 10 × 10 grid of subvoxel MD values (horizontal axes) and subvoxel T₂ values (vertical axes).

The spatially-resolved subvoxel T₁-T₂ and MD-T₂ spectral components from the second TBI case (Subject 3) are shown in Figure 3. Unfiltered (T₁-T₂ and MD-T₂ in Figures 3A, C, respectively) and filtered data (T₁-T₂ and MD-T₂ in Figures 3B, D, respectively) are shown. As before, signal components that were spatially consistent with specific tissue types as a function of T₁, T₂, and MD were revealed.

FIGURE 3

FIGURE 3. Maps of 2D probability density functions (i.e., 2D normalized spectra) from Subject 3 (TBI) of (A) unfiltered and (B) filtered subvoxel T₁-T₂ values reconstructed on a 10 × 10 grid of subvoxel T₁ values (horizontal axes) and subvoxel T₂ values (vertical axes), and maps of (C) unfiltered and (D) filtered subvoxel MD-T₂ values reconstructed on a 10 × 10 grid of subvoxel MD values (horizontal axes) and subvoxel T₂ values (vertical axes).

Figure 4 shows histological images and multidimensional MR-derived injury biomarker maps of the three representative cases. Histological images (red = APP stain) of the control case (Subject 1) show negative APP staining, compared with positive APP staining in the injured samples (Subjects two and 3). We then examine separately the two MRI-derived injury biomarkers, T₁-T₂ and MD-T₂, and show the resulting images obtained using the unfiltered full dataset (as originally published in [22]), the filtered full dataset, and the filtered reduced dataset. In addition, the MRI-derived injury biomarkers obtained by using the unfiltered reduced dataset are shown in Supplementary Figure S1 in the Supplementary Material.

FIGURE 4

FIGURE 4. Histological images and multidimensional MR-derived injury biomarker maps of three representative cases, and under different conditions (left to right: unfiltered, filtered, and filtered and reduced data). Deconvolved histological APP images (co-registered with the MRI) are shown on the left panel, red = APP stain (top to bottom: control, and two TBI cases). All multidimensional injury maps were thresholded at 10% of the maximal intensity and overlaid on grayscale proton density images. Multidimensional injury maps of Subject 1 (control) show absent of significant injury under all experimental conditions. Multidimensional injury maps of Subject 2 (TBI) show substantial injury along the white-gray matter interface under all experimental conditions. Multidimensional injury maps of Subject 3 (TBI) show substantial injury at the bottom of the CC under all experimental conditions.

Visual inspection of the different injury biomarker maps shown in Figure 4 and Supplementary Figure S1 revealed that filtering of the data does not result in loss of the spectral information of interest, and furthermore, the filtered images appear qualitatively of higher quality. Importantly, the data reduction in the case of the filtered data did not significantly affect the resulting injury biomarker maps (Figure 4).

3.2 Evaluation of Performance and Correlation With Histology

Evaluation of filtering performance was based upon the extent of noise reduction and feature preservation, and was quantified by computing the structural similarity index (SSIM) values [33] between the injury biomarker maps under the different experimental conditions (e.g., unfiltered, filtered) and the co-registered APP density histological image as reference. All of the SSIM values are shown in Figure 5. In the context of the current study we are most interested in the ability to accelerate the multidimensional MRI acquisition, and therefore the accuracy and quality of the reduced data cases are of particular importance. Compared with the unfiltered and reduced data injury biomarker maps, the SSIM values of the filtered and reduced data images increased by 11.1, 0.9, and 14.3% for the MD-T₂-based biomarker for Subjects 1 to 3, respectively, and increased by 8.6, 7.7, and 4.6% for the T₁-T₂-based biomarker for Subjects 1 to 3, respectively. All of these increases in SSIM were statistically significant (p < 0.001).

FIGURE 5

FIGURE 5. The structural similarity index (SSIM) values between the injury biomarker images under the different experimental conditions (e.g., T₁-T₂ unfiltered, MD-T₂ filtered reduced data) and the co-registered APP density histological image as reference. The three bars at each condition represent the different Subjects (blue = Subject 1, red = Subject 2, and yellow = Subject 3).

To further evaluate the performance of the NESMA filter and the subsequent data reduction, we performed radiological–pathological correlation analyses with histological APP density and all the investigated MRI parameters under the different experimental conditions (e.g., MD-T₂ unfiltered, T₁-T₂ filtered and reduced). Figure 6 summarizes the association between the investigated MR metrics and the pathological findings in normal WM, cortical GM, and TAI ROIs.

FIGURE 6

FIGURE 6. APP density (% area) from 36 tissue regions (APP-positive regions from each TBI case, WM and GM regions), and its correlation with injury biomarker parameter under different experimental conditions. Individual data points represent the mean ROI value from each post-mortem tissue sample. Scatterplots of the mean (with 95% confidence interval error bars) % area APP and (A) MD-T₂ unfiltered (B) MD-T₂ unfiltered and reduced (C) MD-T₂ filtered (D) MD-T₂ filtered and reduced (E)T₁-T₂ unfiltered (F)T₁-T₂ unfiltered and reduced (G)T₁-T₂ filtered, and (H)T₁-T₂ filtered and reduced, show positive and significant correlation with APP density.

To assess the relationship of the MRI parameters with the degree of injury, all tissue ROIs were grouped together and correlated with the APP density (solid lines in Figure 6). All multidimensional injury biomarker cases under every experimental condition were strongly and significantly positively correlated with the APP density (i.e., % area APP). These correlations illustrate how the multidimensional injury biomarker maps provide “true negative”, in the sense that any region outside of the TAI lesions has zero or close to zero intensity.

Of particular interest was the case of the MD-T₂ injury biomarker, where the unfiltered dataset resulted in relatively scattered correlation (r = 0.565, p < 2 × 10^–4, Figure 6A), which was largely unaffected by the data reduction (r = 0.557, p < 3 × 10^–4, Figure 6B). Marked improvement was observed after filtering the data, which led to significantly tighter correlation between the MRI biomarker and the histological marker (r = 0.812, p < 2 × 10^–9, Figure 6C). As can be expected, the 59.6% reduction in the data amount led to some reduction in the radiological–pathological correlation (r = 0.723, p < 4 × 10^–7, Figure 6D), although still performing better than the unfiltered full MD-T₂ dataset. Furthermore, filtering led to a reduced variance within the ROIs, as seen by the decrease in the 95% confidence intervals (error bars in Figure 6).

The T₁-T₂ injury biomarker derived from the unfiltered dataset had excellent correlation with histological APP density (r = 0.970, p < 1 × 10^–22, Figure 6E), and therefore, it is not surprising that this strong relationship was maintained under all of the investigated experimental conditions (r = 0.963, p < 1 × 10^–20, r = 0.962, p < 1 × 10^–20, and r = 0.969, p < 1 × 10^–21, Figures 6F–H, for unfiltered and reduced, filtered, and filtered and reduced, respectively).

4 Discussion

Here we report the use of the NESMA filter on multidimensional MRI data, in particular voxelwise T₁-T₂ and MD-T₂ spectra in fixed human Corpus Callosum, to remove noise and reduce total scan time. We focused on results from a recent study that identified potential imaging biomarkers of axonal injury pathology from the joint analysis of multidimensional MRI and histopathological data [22]. These axonal injury maps were shown to be significantly and strongly correlated with histological evidence of axonal injury. Reprocessing these data provided an opportunity to test the performance of the NESMA filter and its effect on the accuracy of the injury biomarker maps, relative to the histological reference.

Our findings showed that noise reduction in the multidimensional MRI data using an adaptive nonlocal multispectral filter (i.e., NESMA [20]) improved the accuracy of the resulting injury biomarker maps, and furthermore, allowed for data reduction of 35.7 and 59.6% from the full dataset, which led to using only 36 and 122 images in the T₁-T₂ and MD-T₂ cases, respectively.

Specifically, visual inspection and a side-by-side comparison of the unfiltered and filtered subvoxel T₁-T₂ and MD-T₂ spectral components (Figures 1–3) showed that the filtered maps exhibit lower random variations, in particular at the lower ends of the spectra, and that there was no apparent loss of spectral information. For example, Subject 3 exhibited a relatively focal axonal injury at the bottom of the CC (Figure 4, left panel), captured at the lower end of the T₁-T₂ spectra, which was previously associated with axonal injury [22]. Noticeable noise reduction at these spectral lower ends was observed, which is crucial to the robust identification of axonal injury from these multidimensional MRI data.

Visual inspection of the resulting T₁-T₂ and MD-T₂ injury biomarker maps with respect to co-registered APP histological images suggested improved accuracy after applying the NESMA filter, even after the data was reduced (Figure 4). Quantitative evaluation that compared the SSIM between the co-registered APP histological images and injury biomarker maps derived from unfiltered and filtered reduced datasets showed a significant increase as a result of the filtering across all subjects (Figure 5).

We performed radiological–pathological correlation analyses with histological APP density and all the investigated MRI parameters under the different experimental conditions to assess quantitatively whether and to what extent the proposed approach preserves strong correlations even under substantial data reduction (Figure 6). This analysis indicated that not only the correlations were preserved, but furthermore, they were considerably improved, even after data reduction, as a result of filtering the data. Lastly, our results suggest that the previously proposed [22] adaptive method of locating the injury-associated T₁-T₂-MD spectral signature is robust to noise removal procedures and to data reduction.

Common to all ex vivo human MRI studies, our data include the effects of post-mortem degeneration, fixation and resulting dehydration. Because T₁, T₂, and diffusion dynamics are different in fixed tissue compared with living systems, further investigation will be needed to establish whether and how the axonal injury-related T₁-T₂-MD range of multidimensional magnetic resonance parameters is altered in vivo. In this context it is important to note that our findings are not based on absolute values of T₁, T₂, and MD, which indeed are expected to change in vivo. Instead, all of the injury biomarkers maps are generated using the relative signal fraction of an automatically identified T₁-T₂-MD range that does not depend on the actual values of these parameters [22]. We therefore anticipate that in vivo spectra will be shifted in all T₁, T₂, and MD dimensions compared to our ex vivo findings, however, the distributions of the signal fractions should largely remain similar.

The performance of the NESMA filter does not depend on the particular multidimensional spectra quantification processing pipeline because the filter is applied in the image domain, before its transformed into voxelwise spectra. Here we applied a constrained ℓ₂ regularized inversion framework [13], however, we anticipate that the demonstrated improvement after filtering could be extended to other approaches such as ℓ₁ regularization [34], Monte-Carlo inversion [35, 36], and InSpect [16].

Multidimensional MRI is an emerging approach [37] that is now being applied to address a range of medical conditions such as prediction of pregnancy complications via placenta characterization [9], spinal cord injury [6, 38], prostate cancer [39], breast cancer [40], and axonal injury due to TBI [22]. Recent in vivo proof-of-concept applications of subvoxel T₁-T₂ correlation spectra using 105 images [41] and of subvoxel diffusion-T₁ correlation spectra using 363 [11] and 304 [42] images are promising. Here we showed that accurate and robust subvoxel T₁-T₂ and MD-T₂ correlation spectra can be obtained using only 36 and 122 images, respectively, by using a constrained optimization data processing framework (i.e., MADCO [13]) in conjunction with applying the NESMA filter to reduce noise. A reliable and robust noise removal and consequent acquisition acceleration should further advance the field towards clinically-feasible diagnostic multidimensional MRI protocols.

Data Availability Statement

The datasets generated and analyzed during the current study are available from the corresponding author on reasonable request. Requests to access these datasets should be directed to andlt; basserp@mail.nih.gov.

Ethics Statement

For each case, a next–of-kin or legal representative provided a written consent for donation of the brain for use in research.

Author Contributions

DB: conceptualization, design of the study, methodology, software, investigation, data curation, writing—original draft, writing—review and editing, visualization, supervision, and project administration. MB: conceptualization, design of the study, methodology, software, and writing—review and editing. MK: methodology, investigation, and writing—review and editing. DI: methodology, investigation, and writing—review and editing. DP: investigation, methodology, resources, and writing—review and editing. DLB: design of the study, investigation, resources, and writing—review and editing. PB: conceptualization, methodology, resources, writing—review and editing, and funding acquisition. All authors contributed to the article and approved the submitted version.

Funding

This research was partially supported by a grant from the U.S. Department of Defense, Program Project 308430 Uniformed Services University of the Health Sciences (USUHS). Support for this work also included funding from the U.S. Department of Defense to the Brain Tissue Repository and Neuropathology Program, Center for Neuroscience and Regenerative Medicine (CNRM). DB and MEK were supported by the CNRM Neuroradiology- Neuropathology Correlations Core. MB was supported by the Intramural Research Program of the National Institute on Aging. DI, DPP, and DLB were supported by the CNRM and USUHS. PJB was supported by the Intramural Research Program of the Eunice Kennedy Shriver National Institute of Child Health and Human Development.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank the subjects’ families that consented for brain donations for the better understanding of TBI consequences. The authors thank Patricia Lee, Nichelle Gray and Paul Gegbeh for their valuable technical work. We are grateful to Stacey Gentile, Deona Cooper and Harold Kramer Anderson for their administrative assistance. We thank the TRACK-TBI Investigators (https://tracktbi.ucsf.edu/transforming-research-and-clinical-knowledge-tbi). The opinions expressed herein are those of the authors and are not necessarily representative of those of the Uniformed Services University of the Health Sciences (USUHS), the Department of Defense (DOD), the NIH or any other US government agency.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2021.737374/full#supplementary-material

References

1. English AE, Whittall KP, Joy MLG, Henkelman RM. Quantitative Two-Dimensional Time Correlation Relaxometry. Magn Reson Med (1991) 22:425–34. doi:10.1002/mrm.1910220250

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Hürlimann MD, Flaum M, Venkataramanan L, Flaum C, Freedman R, Hirasaki GJ. Diffusion-relaxation Distribution Functions of Sedimentary Rocks in Different Saturation States. Magn Reson Imaging (2003) 21:305–10. doi:10.1016/s0730-725x(03)00159-0

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Topgaard D. Multidimensional Diffusion MRI. J Magn Reson (2017) 275:98–113. doi:10.1016/j.jmr.2016.12.007

CrossRef Full Text | Google Scholar

4. Benjamini D, Basser PJ. Multidimensional Correlation MRI. NMR Biomed (2020) 33:e4226. doi:10.1002/nbm.4226

PubMed Abstract | CrossRef Full Text | Google Scholar

5. de Almeida Martins JP, Topgaard D. Two-Dimensional Correlation of Isotropic and Directional Diffusion Using NMR. Phys Rev Lett (2016) 116:087601. doi:10.1103/PhysRevLett.116.087601

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Kim D, Doyle EK, Wisnowski JL, Kim JH, Haldar JP. Diffusion‐relaxation Correlation Spectroscopic Imaging: A Multidimensional Approach for Probing Microstructure. Magn Reson Med (2017) 78:2236–49. doi:10.1002/mrm.26629

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Benjamini D, Basser PJ. Water Mobility Spectral Imaging of the Spinal Cord: Parametrization of Model-free Laplace MRI. Magn Reson Imaging (2019) 56:187–93. doi:10.1016/j.mri.2018.12.001

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Benjamini D, Basser PJ. Magnetic Resonance Microdynamic Imaging Reveals Distinct Tissue Microenvironments. NeuroImage (2017) 163:183–96. doi:10.1016/j.neuroimage.2017.09.033

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Slator PJ, Hutter J, Palombo M, Jackson LH, Ho A, Panagiotaki E, et al. Combined Diffusion‐relaxometry MRI to Identify Dysfunction in the Human Placenta. Magn Reson Med (2019) 82:95–106. doi:10.1002/mrm.27733

PubMed Abstract | CrossRef Full Text | Google Scholar

10. de Almeida Martins JP, Tax CMW, Reymbaut A, Szczepankiewicz F, Chamberland M, Jones DK, et al. Computing and Visualising Intra‐voxel Orientation‐specific Relaxation-Diffusion Features in the Human Brain. Hum Brain Mapp (2021) 42:310–28. doi:10.1002/hbm.25224

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Reymbaut A, Critchley J, Durighel G, Sprenger T, Sughrue M, Bryskhe K, et al. Toward Nonparametric Diffusion‐ Characterization of Crossing Fibers in the Human Brain. Magn Reson Med (2021) 85:2815–27. doi:10.1002/mrm.28604

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Benjamini D. Chapter 10. Nonparametric Inversion of Relaxation and Diffusion Correlation Data. In: D Topgaard, editor. Advanced Diffusion Encoding Methods in MRI. Cambridge, UK: Royal Society of Chemistry (2020). p. 278–316. chap. 10. doi:10.1039/9781788019910-00278

CrossRef Full Text | Google Scholar

13. Benjamini D, Basser PJ. Use of Marginal Distributions Constrained Optimization (MADCO) for Accelerated 2D MRI Relaxometry and Diffusometry. J Magn Reson (2016) 271:40–5. doi:10.1016/j.jmr.2016.08.004

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Benjamini D, Basser PJ. Towards Clinically Feasible Relaxation-Diffusion Correlation MRI Using MADCO. Microporous Mesoporous Mater (2018) 269:93–6. doi:10.1016/j.micromeso.2017.02.001

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Bai R, Cloninger A, Czaja W, Basser PJ. Efficient 2D MRI Relaxometry Using Compressed Sensing. J Magn Reson (2015) 255:88–99. doi:10.1016/j.jmr.2015.04.002

CrossRef Full Text | Google Scholar

16. Slator PJ, Hutter J, Marinescu RV, Palombo M, Jackson LH, Ho A, et al. Data-Driven Multi-Contrast Spectral Microstructure Imaging with InSpect: INtegrated SPECTral Component Estimation and Mapping. Med Image Anal (2021) 71:102045. doi:10.1016/j.media.2021.102045

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Hutter J, Slator PJ, Christiaens D, Teixeira RPAG, Roberts T, Jackson L, et al. Integrated and Efficient Diffusion-Relaxometry Using ZEBRA. Sci Rep (2018) 8:15138. doi:10.1038/s41598-018-33463-2

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Manhard MK, Stockmann J, Liao C, Park D, Han S, Fair M, et al. A Multi‐inversion Multi‐echo Spin and Gradient echo echo Planar Imaging Sequence with Low Image Distortion for Rapid Quantitative Parameter Mapping and Synthetic Image Contrasts. Magn Reson Med (2021) 86:866–80. doi:10.1002/mrm.28761

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Bouhrara M, Bonny J-M, Ashinsky BG, Maring MC, Spencer RG. Noise Estimation and Reduction in Magnetic Resonance Imaging Using a New Multispectral Nonlocal Maximum-Likelihood Filter. IEEE Trans Med Imaging (2017) 36:181–93. doi:10.1109/tmi.2016.2601243

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Bouhrara M, Reiter DA, Maring MC, Bonny J-M, Spencer RG. Use of the NESMA Filter to Improve Myelin Water Fraction Mapping with Brain MRI. J Neuroimaging (2018) 28:640–9. doi:10.1111/jon.12537

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Bouhrara M, Lee DY, Rejimon AC, Bergeron CM, Spencer RG. Spatially Adaptive Unsupervised Multispectral Nonlocal Filtering for Improved Cerebral Blood Flow Mapping Using Arterial Spin Labeling Magnetic Resonance Imaging. J Neurosci Methods (2018) 309:121–31. doi:10.1016/j.jneumeth.2018.08.018

CrossRef Full Text | Google Scholar

22. Benjamini D, Iacono D, Komlosh ME, Perl DP, Brody DL, Basser PJ. Diffuse Axonal Injury Has a Characteristic Multidimensional MRI Signature in the Human Brain. Brain (2021) 144:800–16. doi:10.1093/brain/awaa447

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Johnson VE, Stewart W, Smith DH. Axonal Pathology in Traumatic Brain Injury. Exp Neurol (2013) 246:35–43. doi:10.1016/j.expneurol.2012.01.013

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Pas K, Komlosh ME, Perl DP, Basser PJ, Benjamini D. Retaining Information from Multidimensional Correlation MRI Using a Spectral Regions of Interest Generator. Sci Rep (2020) 10:3246. doi:10.1038/s41598-020-60092-5

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Provencher SW. A Constrained Regularization Method for Inverting Data Represented by Linear Algebraic or Integral Equations. Comp Phys Commun (1982) 27:213–27. doi:10.1016/0010-4655(82)90173-4

CrossRef Full Text | Google Scholar

26. Kroeker RM, Henkelman MR. Analysis of Biological NMR Relaxation Data with Continuous Distributions of Relaxation Times. J Magn Reson (1969) 69:218–35. doi:10.1016/0022-2364(86)90074-0

CrossRef Full Text | Google Scholar

27. Song Y-Q, Venkataramanan L, Hürlimann MD, Flaum M, Frulla P, Straley C. T1-T2 Correlation Spectra Obtained Using a Fast Two-Dimensional Laplace Inversion. J Magn Reson (2002) 154:261–8. doi:10.1006/jmre.2001.2474

CrossRef Full Text | Google Scholar

28. Mitchell J, Chandrasekera TC, Gladden LF. Numerical Estimation of Relaxation and Diffusion Distributions in Two Dimensions. Prog Nucl Magn Reson Spectrosc (2012) 62:34–50. doi:10.1016/j.pnmrs.2011.07.002

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Celik H, Bouhrara M, Reiter DA, Fishbein KW, Spencer RG. Stabilization of the Inverse Laplace Transform of Multiexponential Decay through Introduction of a Second Dimension. J Magn Reson (2013) 236:134–9. doi:10.1016/j.jmr.2013.07.008

CrossRef Full Text | Google Scholar

30. Labadie C, Lee JH, Vetek G, Springer CS. Relaxographic Imaging. J Magn Reson Ser B (1994) 105:99–112. doi:10.1006/jmrb.1994.1109

CrossRef Full Text | Google Scholar

31. Macenko M, Niethammer M, Marron JS, Borland D, Xiaojun Guan JT. A Method for Normalizing Histology Slides for Quantitative Analysis. In IEEE International Symposium on Biomedical Imaging. From Nano to MacroIEEE (2009) p. 1107–10. doi:10.1109/ISBI.2009.5193250

CrossRef Full Text | Google Scholar

32. Ruifrok AC, Johnston DA. Quantification of Histochemical Staining by Color Deconvolution. Anal Quant Cytol Histol (2001) 23:291–9.

PubMed Abstract | Google Scholar

33. Wang Z, Bovik AC, Sheikh HR, Simoncelli EP. Image Quality Assessment: From Error Visibility to Structural Similarity. IEEE Trans Image Process (2004) 13:600–12. doi:10.1109/tip.2003.819861

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Reci A, Sederman AJ, Gladden LF. Obtaining Sparse Distributions in 2D Inverse Problems. J Magn Reson (2017) 281:188–98. doi:10.1016/j.jmr.2017.05.010

CrossRef Full Text | Google Scholar

35. Prange M, Song Y-Q. Quantifying Uncertainty in NMR Spectra Using Monte Carlo Inversion. J Magn Reson (2009) 196:54–60. doi:10.1016/j.jmr.2008.10.008

PubMed Abstract | CrossRef Full Text | Google Scholar

36. de Almeida Martins JP, Topgaard D. Multidimensional Correlation of Nuclear Relaxation Rates and Diffusion Tensors for Model-free Investigations of Heterogeneous Anisotropic Porous Materials. Sci Rep (2018) 8:2488. doi:10.1038/s41598-018-19826-9

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Slator PJ, Palombo M, Miller KL, Westin C, Laun F, Kim D, et al. Combined Diffusion-Relaxometry Microstructure Imaging: Current Status and Future Prospects. Magn Reson Med (2021). doi:10.1002/mrm.28963

CrossRef Full Text | Google Scholar

38. Benjamini D, Hutchinson EB, Komlosh ME, Comrie CJ, Schwerin SC, Zhang G, et al. Direct and Specific Assessment of Axonal Injury and Spinal Cord Microenvironments Using Diffusion Correlation Imaging. NeuroImage (2020) 221:117195. doi:10.1016/j.neuroimage.2020.117195

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Zhang Z, Wu HH, Priester A, Magyar C, Afshari Mirak S, Shakeri S, et al. Prostate Microstructure in Prostate Cancer Using 3-T MRI with Diffusion-Relaxation Correlation Spectrum Imaging: Validation with Whole-Mount Digital Histopathology. Radiology (2020) 296:348–55. doi:10.1148/radiol.2020192330

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Naranjo ID, Reymbaut A, Brynolfsson P, Lo Gullo R, Bryskhe K, Topgaard D, et al. Multidimensional Diffusion Magnetic Resonance Imaging for Characterization of Tissue Microstructure in Breast Cancer Patients: A Prospective Pilot Study. Cancers (2021) 13:1606. doi:10.3390/cancers13071606

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Kim D, Wisnowski JL, Nguyen CT, Haldar JP. Multidimensional Correlation Spectroscopic Imaging of Exponential Decays: From Theoretical Principles to In Vivo Human Applications. NMR Biomed (2020) 33. doi:10.1002/nbm.4244

CrossRef Full Text | Google Scholar

42. Avram AV, Sarlls JE, Basser PJ. Whole-Brain Imaging of Subvoxel T1-Diffusion Correlation Spectra in Human Subjects. Front Neurosci (2021) 15:671465. doi:10.3389/fnins.2021.671465

PubMed Abstract | CrossRef Full Text | Google Scholar