Mahalanobis distance tractometry (MaD-Tract) - a framework for personalized white matter anomaly detection applied to TBI.

Imaging-based quantitative measures from diffusion-weighted MRI (dMRI) offer the ability to non-invasively extract microscopic information from human brain tissues. Group-level comparisons of such measures represent an important approach to investigate abnormal brain conditions. These types of analyses are especially useful when the regions of abnormality spatially coincide across subjects. When this is not true, approaches for individualized analyses are necessary. Here we present a framework for single-subject multidimensional analysis based on the Mahalanobis distance. This is conducted along specific white matter pathways represented by tractography-derived streamline bundles. A definition for abnormality was constructed from Wilk's criterion, which accounts for normative sample size, number of features used in the Mahalanobis distance, and multiple comparisons. One example of a condition exhibiting high heterogeneity across subjects is traumatic brain injury (TBI). Using the Mahalanobis distance computed from the three eigenvalues of the diffusion tensor along the cingulum, uncinate, and parcellated corpus callosum tractograms, 8 severe TBI patients were individually compared to a normative sample of 49 healthy controls. For all TBI patients, the analyses showed statistically significant deviations from the normative data at one or multiple locations along the analyzed bundles. The detected anomalies were widespread across the analyzed tracts, consistent with the expected heterogeneity that is hallmark of TBI. Each of the controls subjects was also compared to the remaining 48 subjects in the control group in a leave-one-out fashion. Only two segments were identified as abnormal out of the entire analysis in the control group, thus the method also demonstrated good specificity.

Introduction

Capturing individual variability of quantitative neuroimaging measures in meaningful ways is relevant for gaining insight into the pathophysiology of highly heterogeneous neurological conditions such as traumatic brain injury (TBI) and autism. One important example of single-subject analysis frameworks can be found in the work published by Kim et al. (2013) in the study of TBI. The authors presented a voxel-based analysis (VBA) method based on estimating Z-scores between the fractional anisotropy (FA) of a TBI patient and a healthy control group for each voxel in the brain. The method was able to find unique patterns of voxels with abnormally high or low FA in TBI patients. However, the method displayed specificity limitations since a considerable number of abnormal voxels were also flagged in healthy controls.

Another study by Dean et al. (2017) found a dramatic improvement in the ability to distinguish autistic individuals from a typically developing reference group when using a multivariate approach. In their work, multiple DTI parameters extracted as averages from white matter regions were used to compute the Mahalanobis distance between an individual and the typically developing group. The Mahalanobis distance (Mahalanobis, 1936) has often been used as a measure of class dissimilarity and as a tool for anomaly detection in multivariate normal data in a wide array of applications including the study of autism (Dean et al., 2017), TBI (Taylor et al., 2020), epilepsy (Gyebnár et al., 2019; Hong et al., 2014, 2017; Liu et al., 2016; Morgan et al., 2021; Owen et al., 2020; Pressl et al., 2019; Weng et al., 2020), and others (Jin et al., 2012; Lin et al., 2010; Patil et al., 2015; Wang et al., 2013; Wu and Zhang, 2006; Zhang et al., 2016).

We previously introduced a voxel-wise implementation of the Mahalanobis distance (MaD-Vox) using DTI parameters to compare an individual to a reference group (Guerrero et al., 2018). This study was based on simulated data and demonstrated the ability to capture individual variability in DTI parameters at the voxel level with high specificity. However, when applied to real human brain data, the voxel-based Mahalanobis method was highly susceptible to inter-subject image misalignment. Also, multiple-comparisons correction greatly reduced the statistical power of this method as with most VBA approaches. A recently published description of a VBA Mahalanobis distance method in a study of epilepsy reported similar limitations (Gyebnár et al., 2019).

The work described here outlines a novel tractometry-based (Yeatman et al., 2012) multivariate computational and statistical framework for testing whether an individual deviates significantly from a reference group at the white matter tract level. The method uses the Mahalanobis distance as a multivariate metric for abnormality testing. The testing is performed at individual segments along specific white matter pathways. This reduces the number of performed tests by orders of magnitude compared to VBA but retains better spatial specificity than the ROI-averaging approach. Also, the requirements for spatial image alignment should be less stringent than required for VBA.

Here, the new MaD-Tract framework was demonstrated in a study of white matter microstructure in a cohort of children with severe traumatic brain injuries (Fig. 1). TBI occurs when an external force acts on the head or body leading to neuropathologic damage and functional impairment. TBI may cause life-long health deficits or death, and it represents a major public health concern involving serious socioeconomic burdens nationally as well as worldwide (James et al., 2019; Langlois et al., 2006; Thurman et al., 1999). TBI is the most common form of acquired brain injury and prevalent across pediatric, adult, and aging populations with leading causes being military blasts, vehicular accidents, falls, assaults, and sport activities (Langlois et al., 2006). Clinical manifestation of TBI is extremely heterogeneous depending on mechanism, location, severity, existing comorbidities, and time since injury.

Fig. 1.

(a) Acute anatomical brain scans for 5 of the 8 patients included in this study reveal severe traumatic brain injuries. Scanning during the acute phase following injury was not available for the three additional patients (A, B, and R). (b) Follow-up FLAIR images for all 8 patients in this study acquired 1–2 years post injury.

This study evaluated the MaD-Tract framework on a subset of white mater tracts, specifically the cingulum bundles, uncinate fasciculi, and corpus callosum. These pathways are critical parts of networks that are known to be implicated in the clinical and behavioral outcomes following TBI (Anderson and Catroppa, 2007; Ewing-Cobbs et al., 2008; Fay et al., 2009; Hillary et al., 2011; Johnson et al., 2011; Juranek et al., 2012; Hillary et al., 2011;Anderson and Catroppa, 2007; Johnson et al.,2011; Juranek et al., 2012; Fay et al., 2009; Ewing-Cobbs et al., 2008). In particular, the cingulum bundles include important connections of the default mode network and hippocampal memory network. The uncinate fasciculus is an important pathway in the fronto-limbic network. The corpus callosum is responsible for interhemispheric communication and plays an important role in overall brain integration.

Methods

Participants

Imaging research participants were 11–18 years of age and enrolled in the Approaches and Decisions in Pediatric TBI (ADAPT) trial, which included roughly 1000 children (from birth to 18 years) with severe TBI defined as a post-resuscitation Glasgow Coma Scale (GCS) of 8 or less (Ferrazzano et al., 2019; Sarnaik et al., 2018). A group of twenty-two ADAPT participants were recruited and enrolled into this MRI study ~1–2 years post injury. To best account for site-to-site variation in diffusion measurements, only eight participants from six sites with diffusion calibration phantom scans were included (4 males, 4 females) between the ages of 11.6 and 18.9 years (Mean 15.7 ± 2.1 years). Acute clinical scans for five out of the eight patients are shown in Fig. 1a and follow-up FLAIR images for all eight patients are shown in the Fig. 1b. A typically developing control cohort without history of TBI or neuropsychiatric diagnoses was also recruited for imaging and neurocognitive testing at the University of Wisconsin – Madison (UW). The control group consisted of 49 subjects (24 males, 25 females) between the ages of 9.0 and 18.0 years (mean 13.45 ± 2.8) at the time of MRI scanning. The study was approved by the institutional review board at the University of Wisconsin – Madison and all participating sites. Informed consent was obtained from the subject or legal guardian when appropriate.

Brain imaging

Brain imaging was performed for each participant using 3T MRI standardized neuroimaging protocols across the sites. T1-weighted (T1w), T2-weighted (T2w), T2-weighted FLAIR, T2* - weighted, diffusion tensor, and resting state functional images were obtained for each subject. Manufacturer-specific protocols were emulated after protocols used in the multi-site Transforming Research and Clinical Knowledge in TBI (TRACK-TBI) study.

Prior to subject enrollment, each site was provided with the scanner-specific protocol to be implemented on their system. A scanning procedure manual was developed for the study and disseminated to all participating sites. Prior to human data collection, sites were required to collect phantom data using the provided protocol, which were then sent to UW to verify protocol compliance. Once the site imaging protocol was approved, sites enrolled adolescent TBI participants for scanning. Imaging was performed 12–25.5 months post injury with a mean interval between injury and MRI scanning of 20 ± 4.44 months.

For this analysis, scans that were used included T1w, T2w, and DTI. Structural T1w imaging was performed using a 3D inversion-recovery prepared sequence with a rapid gradient echo readout (i.e., MP-RAGE on Siemens and Philips; BRAVO on GE). The protocol prescription was 3D sagittal images with a 256 × 256 matrix over a 256 mm field of view and 192 slices that were 1 mm thick (1 mm isotropic resolution), frequency encoding in Superior/Inferior direction, in-plane parallel imaging with an acceleration factor of 2. On GE scanners, the inversion time, TI, was set to 450 ms. For Siemens and Philips scanners the TI was 900 ms.

Structural T2w imaging was performed using a 3D fast spin-echo sequence (i.e., SPACE on Siemens; VISTA on Philips; CUBE on GE). The protocol prescription was 3D sagittal images with 2 averages, a 256 × 256 matrix over a 256 mm field of view and 192 slices that were 1 mm thick (1 mm isotropic resolution), frequency encoding in Superior/Inferior direction, in-plane parallel imaging with an acceleration factor of 2. TR/TE was 2500/95 ms for GE, 2500/256 ms Phillips, 2500/398 ms for Siemens.

Diffusion tensor imaging was performed with a single-shot spin-echo echo-planar imaging (EPI) pulse sequence. The protocol consisted of 2D sagittal images, a 96 × 96 matrix, 240 mm FOV, 64 slices, 2.5 mm isotropic resolution, phase encoding in Anterior/Posterior direction, and parallel acquisition with a geometric reduction factor of 2. Diffusion set-tings included diffusion encoding along 64 non-collinear directions with b-value of 1300s/mm2, 8 vol with no diffusion weighting (b = 0 s/mm2). TR/TE was 8500/82 ms for GE, 8500/103 ms Phillips, 9000/82 ms for Siemens.

Data processing

Distortion, translation and rotation from bulk head motion and eddy currents were corrected using the eddy tool (Andersson and Sotiropoulos, 2016) in FSL (version 6.0) with outlier replacement enabled (Andersson et al., 2016). Gradient directions were corrected for image rotations from image registration (Leemans and Jones, 2009). Since field maps were not available, EPI-related geometric distortions were corrected using an in-house image spatial normalization approach, which aligned the diffusion scan images to the structural T2- and T1-weighted images using image co-registration constrained along the phase encoding direction of the diffusion images. The constrained registration in the phase encoding direction was enabled by the use of ‘antsRegistration’ in ANTs (Avants et al., 2011).

Diffusion tensors were estimated for each voxel using the robust estimation of tensors by outlier rejection (RESTORE) algorithm as part of the diffusion imaging in python (DIPY) software package (Garyfallidis et al., 2014). Eigenvalue maps (λ1, λ2, λ3) were generated from the voxel-wise estimates of the diffusion tensor.

Fiber Orientation Distribution Function (fODF) maps were estimated for individual subjects using the MRtrix3 software package (Tournier et al., 2019). White matter response functions for all control subjects were first estimated using dwi2response with the multi-tissue ‘dhollander’ algorithm (Dhollander et al., 2016). These were then averaged across subjects in order to estimate a mean white matter response function. The average white matter response function was then utilized for estimating the fODF maps by constrained spherical deconvolution (CSD) for both control and TBI scans using dwi2fod with the ‘csd’ algorithm (Tournier et al., 2007).

Prior to segmentation of white matter streamline bundles, multi-subject spatial normalization of the fODF maps was performed using diffeomorphic deformations (Raffelt et al., 2011, 2012) to a study-specific fODF template. The template was constructed from all the subjects in the control group using the ‘population_template’ command in Mrtrix3 (Tournier et al., 2019). The subject-to-template transformations were also applied to the DTI scalar maps.

Delineation of white mater fiber bundles was carried out in template space using TractSeg (Wasserthal et al., 2018), which uses a pre-trained convolutional neural network to create region-specific tractograms from fODF peaks. The fODF peaks were created using the sh2peaks tool from MRTrix3 (Jeurissen et al., 2013; Tournier et al., 2019), using the spatially normalized fODF maps as input and a maximum of three peaks per voxel. TractSeg was then used to segment fiber-bundle regions as binary masks as well as bundle start- and end-regions. Tract orientation maps (TOMs) (Wasserthal et al., 2018) were also segmented for each fiber-bundle region from the whole-brain peak maps. Probabilistic fiber tracking was then performed within each TOM using the segmented masks for seeding and the start- and end- regions as inclusion points in order to create bundle specific tractograms.

All reconstructed tracts were inspected visually for defects. However, upon inspection, it was found that TBI subjects S and D with lesions affecting large portions of the left hemisphere did not have any streamlines for the left uncinate bundles. Patient S also had failed reconstructions for CC_1, CC_3, CC_4. Additionally, Patient D had failed reconstruction of CC_4. As a result, these tracts were not included in the analyses.

Finally, tractometry (Yeatman et al., 2012) was conducted on the de-lineated tracts using the algorithm described by Chandio et al. (2019) to generate DTI parameter tract profiles. Tract profiles are represented by vector of mean DTI parameter values sampled from 20 equally distanced segments across the tract, excluding the distal ends of the tracts as these are more likely to be poorly registered across subjects (Yeatman et al., 2012).

The mahalanobis distance

For a multivariate normal random vector = (x1, x2, x3, … , x) of dimension p (also referred to as the feature vector) the multivariate normal probability density function is given by where is the means vector of the distribution, Σ is the distribution covariance matrix, and T denotes vector transpose. The sample means vector and covariance matrix for a single group of individuals may be estimated, respectively, as and where ^ represents an estimated value and N is the number of samples. As a general rule-of-thumb, reliable estimations of the inverse covariance matrix require at least 10 observations per dimension of the feature vector. After substituting the sample means vector and covariance matrix estimates, Eq. (1) becomes

The squared Mahalanobis distance, (MaD)2, is equal to

MaD encodes information about the separation of a measurement from the population means relative to the spread of the distribution about the mean. MaD values reflect the degree that a measurement is deviant relative to a defined normative distribution. An estimate for MaD can be obtained from a population sample by using estimates of the means vector and covariance matrix as

The distribution of (MaD)2 may be approximated by the Chi-squared distribution, but is better approximated by the F distribution for small samples sizes (Penny, 1996). The critical value of MaD for detection of abnormalities may be defined using Wilk’s criterion (Penny, 1996; Wilks, 1963) as which uses an F statistic at a Bonferroni corrected α and degrees of freedom determined by reference sample size N and size of feature vector p, and makes it possible to correct for N and p.

MaD-Tract

In the MaD-Tract ([M]halanobis [D]istance [Tract]ometry) framework, MaD is computed along a specific white matter pathway as depicted in the flow chart in Fig. 2. The first step of this framework is to conduct DTI tractometry (Yeatman et al., 2012), which, as previously discussed, consists of mapping DTI parameters along specific streamline bundles, partitioning the bundle into several segments guided by a centroid streamline, and extracting average parameter values from each segment in order to create a vector of segment-averages. In the MaD-Tract framework, multiple quantitative imaging measures are profiled along a tract for a reference group as well as for the individual to be analyzed. The reference group measures are used to define the normative mean and covariance for each tract partition of the tract. For the subject being evaluated, the MaD is calculated using Eq. (6) at each tract partition. A tract partition is considered abnormal if MaD exceeds the critical value in Eq. (7).

Fig. 2.

Steps of the MaD-Tract framework. A. Tractometry: a streamline bundle is partitioned into a number of segments guided by a centroid streamline; DTI eigenvalue maps are also mapped onto the streamline bundle; the partitioned bundle and parameters are combined to create parameter profiles. B. Data from the profiles is structured into a reference group and a single subject profiles. C. At each individual partition (or segment) the Mahalanobis distance between the individual and the group is computed. A Mahalanobis distance profile along the tract is produced.

MaD-Tract was used to evaluate each of the selected tractograms for the 8 severe TBI patients in comparison to the reference control group of 49 subjects. Additionally, in a leave-one-out fashion, each of the controls was compared to the rest of the group by removing it before estimating the mean vector and covariance matrix. The principal diffusivities (eigenvalues: λ1, λ2, λ3) of the diffusion tensors were used for Mahalanobis distance evaluations. Data structuring and estimation of the Mahalanobis distance (including means vectors and inverse covariance matrices) was conducted with in-house developed scripts in the MAT-LAB computer language (Mathworks, Inc). Here testing was performed at a significance level of 0.05 (Bonferroni corrected α = 3.7 × 10−6 for 8 TBI patients, 49 controls, 11 streamline bundles with 18 segments per bundle, n = 49, and p = 3 (the number of features for 3 eigenvalues), resulting in a critical MaD > 6.38.

Multi-site data harmonization

In order to account for potential systematic effects related to site on the estimated tensor eigenvalues, we collected diffusion phantom scans at the 6 sites. The scanned phantom was the Diffusion Phantom Model 128 (High Precision Devices, Inc, Boulder, CO) developed by the National Institute of Standards and Technology (NIST) and the Radiological Society of North America (RSNA)’s Quantitative Imaging Biomarker Alliance (QIBA). As shown in Fig. 3, the phantom contains 13 vials, two for each of 5 concentrations (10, 20, 30, 40, 50%) of polyvinylpyrrolidone (PVP) in an aqueous solution for modulating the isotropic diffusivity of water protons, plus three vials with no PVP (i.e. 0% concentration) (Boss et al., 2015). One of the 0% PVP vials is at the center of the housing and one at each of two concentric rings.

Fig. 3.

Harmonization relationships. These plots show linear regressions estimated between measurements of diffusion coefficient for the 5 PVP concentrations between Site I and each of the other 5 sites in which the phantom was scanned. Also shown is a picture of the NIST PVP diffusion phantom with varying PVP concentration vials. Note there are two vials per concentration, each placed in one of two concentric circles (except for the 0%, which has 3 vials, one in the center and one in each of the circles).

The phantom was scanned at 0 °C, which enabled direct comparison of estimated diffusion measurements to standardized diffusion coefficient values and also between sites. This wa achieved by placing the phantom in an ice-water bath for a minimum of 2.5 h. Temperature recordings from each site of the phantom before and after scanning were recorded. Preparation instructions were developed at UW-Madison and distributed to the other sites. Phantom measurements were obtained from each scanner vendor (GE, Phillips, Siemens). The human study DTI protocol was used for the phantom studies and data was processed the same as the human brain scans. Regions defined as cross sections near the center of each vial were used to extract average diffusion values.

The diffusion phantom data collected across sites was used to estimate harmonization relationships with respect to the control group site. These are shown in Fig. 3. These site-to-site linear relationships were derived from using the multiple PVP concentration vials and their estimated tensor diagonal diffusivities and fitting a linear model that predicts the diffusivity of a given site as a function of the diffusivity estimated in the reference site. Given that these are isotropic media and all three tensor diagonal elements (e.g. Dxx, Dyy, Dzz) appear to follow similar trends, a single linear regression was estimated for each site. These relationships demonstrate remarkable consistency of the diffusion measurements across sites. Nonetheless, before computing the Mahalanobis distance, these linear regression relationships were used to harmonize the eigenvalue parameter profiles of each site and the control group.

Results

The MaD-Tract results for the two patients scanned at the same site as the control subjects are presented in Figs. 4-6. More specifically, Fig. 4 displays for patient A all the 8 skeletonized and partitioned tracts with detected anomalies shown in red. The white matter tracts where anomalies were detected are CC_4 and CC_5. The profiles along these two tracts for Mahalanobis distance as well as for the eigen values are also shown. These plots contain the profiles of both the patient and the control subjects. For this specific case, the profiles shown in these plots are helpful in conveying the potential for the multivariate approach to increase the ability to detect deviations from normal. The univariate profiles do not seem to provide enough information by themselves for classifying certain segments as abnormal. On the other hand, the Mahalanobis distance profiles show some of tract partitions are clearly deviant. These are the segments that exceed the threshold calculated using Eq. (7), shown in the profile plots as the shaded rectangular area.

Fig. 4.

Patient-specific microstructural anomaly detection with MaD-Tract. These results show skeletonized tracts for patient A in the top panel with abnormalities shown in red. The segments flagged as abnormal exceed a predetermined abnormality index depicted as the shaded gray area in the MaD plots (lower panel). The distribution of MaD as well as the univariate parameter profiles for the reference group are also shown. For reference, the skeletonized tracts color coded by direction embedded in the T1w volume are shown in the top left and right corners.

Fig. 6.

Anomalous tract segments on T1w for patient B. The segments flagged as significantly deviant in several tracts for patient B are shown here overlaid or embedded in the subject T1w volume. Visibly atrophied areas (yellow arrows) are observable in the anatomical scan obtained approximately 1year after the sustained severe brain injuries. Some of the detected abnormalities are in the vicinity of these regions. However, anomalies also occur in the absence atrophy observable on the T1w contrast(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.).

In Fig. 5 the abnormal tract partitions are displayed on the T1w anatomical image. Visualization of the anomalies on anatomy facilitates further examination of where MaD-tract detected abnormalities are found with respect to, the location of the original injury. It is worth noting that these are the partitions of the centroid or skeletonized tract and only represent a general trajectory of the tract. Thus, the location of a given MaD-Tract anomaly with respect to a specific anatomical feature is approximate. This is similar to how tract-based spatial statistics (Smith et al., 2006) are typically displayed, except here the skeleton is specific to a white matter tract instead of the whole brain.

Fig. 5.

Anomalous tract segments on T1w. The segments flagged as significantly deviant in CC_4 and CC_5 for patient A are shown here overlaid on the subject T1w volume. Though ventricles are abnormally enlarged, no visibly atrophied areas are observable in the anatomical scan obtained approximately 1year after the sustained severe brain injuries.

The results for patient B are summarized in Fig. 6. This shows only the tracts with detected abnormalities: CC_1, CC_3, CC_5, CC_6, CC_7, and right UF. Compared to patient A, the number of detected abnormal tract partitions is larger and more widespread for patient B. Overlaying the abnormal segments on the T1w image shows some of them occur near areas of atrophied tissue such as in the right uncinate as well as in right hemispheric deep white matter near CC_3. However, some anomalies occur in the absence of visibly abnormal tissue intensities on the T1w image. For example, the segment identified as abnormal in CC_1 of the corpus callosum.

The results corresponding to the multi-site component patients are displayed in Figs. 7 and 8. More specifically, these figures show the Mahalanobis distance parameter profiles for patients and controls. For patients having tract partitions identified as abnormal, the entire profile is highlighted in order to differentiate from those whose all tract partitions are normal. For example, in CC_2 patients C, D, S, and V were found to have abnormalities, but patients R and M did not have any detected. For reference, the Mahalanobis distance for each control subject as well as the group average are also plotted. Figures for these patients and their detected anomalies as overlays on T1w images are provided in the Supplement. The group level display of results in Figs. 7, 8 as well as the figures in the Supplement, reveal the individual variability of abnormalities captured by MaD-Tract. The occurrence of anomalies is both widespread and mostly inconsistent across tracts and participants.

Fig. 7.

Plots of MaD profiles for the cingulum and uncinate bundles. The profiles of tracts with abnormalities are highlighted in subject-specific colors, those who do not present abnormalities are represented by dots (black: control; subject-specific color: TBI) in the shaded region. The shaded region represents the region of normalcy, where the upper bound is set by the critical MaD value of 6.38. Subject-specific coloring of profiles with abnormality facilitates the traceability of the profiles.

Fig. 8.

Plots of MaD profiles for the corpus callosum bundles. The profiles of tracts with abnormalities are highlighted in subject-specific colors, those who do not present abnormalities are represented by dots (black: control; subject-specific color: TBI) in the shaded region. The shaded region represents the region of normalcy, where the upper bound is set by the critical MaD value of 6.38. Subject-specific coloring of profiles with abnormality facilitates the traceability of the profiles.

Finally, a summary of the results for all patients is included in Table 1. While all tracts were flagged, Table 1 shows that the corpus callosum as a whole had the highest frequency of detected anomalies.

Table 1

Summary of results. TBI patients sorted by Site, scanner Vendor, and scanner Model. For each white matter bundle, a mark (‘ x’) is given if at least one segment along the tract was flagged as abnormal.

Site	Vendor	Patient	CORPUS CALLOSUM							CINGULUM		UNCINATE		Number ofbundles withabnormalitiesper patient
			CC_1	CC_2	CC_3	CC_4	CC_5	CC_6	CC_7	left	right	left	right
I	GE	A				×	×							2
I	GE	B	×		×		×	×	×				×	6
II	GE	C	×	×	×		×	×	×	×	×	×	×	10
III	Philips	D	×	×	×	×	×	×	×	×				8
IX	Siemens	M				×	×	×		×				4
XI	Siemens	R	×											1
XI	Siemens	S		×			×	×	×	×	×		×	7
XIII	Siemens	V	×	×		×	×		×	×	×	×	×	9
Number of patients with abnormalities per bundle			5	4	3	4	7	5	5	5	3	2	4

Discussion

In this work we present a multivariate framework for single subject analysis of dMRI measures within a specific white matter tract. The computational and statistical approach captures individual variability by estimating the Mahalanobis distance of multiple variables between a subject and a reference group. Abnormal Mahalanobis distance values are determined according to Wilk’s criterion, which enables the ability to correct for number of dimensions as well as reference group size. Multiple comparisons are addressed by Bonferroni correction of the desired significance level and feeding the result to Wilk’s test for computing a critical Mahalanobis distance value. We refer to the framework as MaD-Tract. The Mahalanobis distance as estimated here exploits existing correlations between the tensor diffusivities through the estimation of their covariance matrix, which may work as an additional source of information and may in turn lead to improved detection of alterations in the diffusional processes measured with DTI. Additionally, the Mahalanobis distance may incorporate other microstructure metrics derived from non-diffusion data such as relaxometry parameters (e.g. R1, R2) or other modalities such as quantitative PET.

MaD-Tract in severe TBI

Implementation of MaD-tract on eight severe TBI cases highlighted individual variability of detected abnormalities. The investigated tracts included the cingulum, uncinate, and corpus callosum, which are expected to be implicated in chronic neurologic deficits in severe TBI. While MaD-Tract detected abnormalities in all cases, their number and location varied across subjects and were more widespread in some cases compared to others. These findings reflect the expected heterogeneity of TBI and highlight the value of single-subject investigation strategies. Nonetheless, in this small severe TBI group the most affected tract in terms of number of detected abnormalities was the corpus callosum, which has been shown to result in global neurologic disfunction when damaged. This finding is consistent with other TBI studies that have also found alterations in FA and MD in the corpus callosum (Rutgers et al., 2008; Wang et al., 2008). Additionally, findings that the corpus callosum undergoes greater shear forces after trauma have been reported in a computer model of fall-induced TBI (Ghajari et al., 2017). Moreover, chronic moderate to severe TBI sequelae have been related to atrophy in the corpus callosum (Kraus et al., 2007). Alterations in FA, MD, or the Mahalanobis distance as estimated here, point to changes in the individual diffusivity values of the diffusion tensor. Though the pathophysiology of TBI at the neuronal level is complex, a growing body of studies have demonstrated that changes at the cellular level induce changes in the magnitude of water diffusion detectable with dMRI. For example, changes in radial diffusivities are thought to be related to changes in myelination, while increases in axial diffusivity are hypothesized to reflect pathology in the axon itself (Kraus et al., 2007). Both types of alterations are involved in diffuse axonal injury (Maxwell et al., 1999; Smith et al., 2003). MaD itself is more agnostic to mechanistic changes and is potentially sensitivity to changes in any of the eigenvalues or combinations of those values.

The profound injuries in the severe TBI study almost guarantees that differences will be found when compared with the healthy group. In some cases, anatomical images such as T1w are able to show focal lesions, particularly in severe cases with gross visible injuries, and the univariate tract profiles in those cases show obvious deviation from the control group. With MaD-tract, abnormalities were found both near and far from visibly atrophied tissue, indicating the method has the ability to detect not only obvious damage but also deviations in regions that may otherwise appear normal. Moreover, when comparing each of the control participants to the other controls in the leave-one-out analysis, only two individual tract sections were flagged as abnormal, demonstrating high specificity of MaD-Tract to severe TBI. It is envisioned and expected that this multivariate framework will prove useful in detecting abnormalities that are more subtle in conditions such as in mild TBI autism or Alzheimer’s disease.

Normative group modeling

The approach described in this work bears similarities to new normative modeling approaches as those presented in Marquand et al. (2016). Such methods model variation either across a study population that includes both healthy and pathological cohorts or across a single large healthy group. In Marquand et al. (2016), Gaussian process regression is utilized to estimate a normative model over a large population sample that links biological variables (as those derived from neuroimaging) to clinically relevant measures (such as trait scores). Along the trajectory defined by the model, confidence centiles are estimated for quantifying ranges of goodness-of-fit between an observation and the normative model. For each observation in the normative population, variation from normal is captured by a Z-score calculated using the estimated centiles. Given the model and its characterized variation, a new observation can be fit to the model and its Z-score with respect to the normal trajectory can be calculated. Extreme value statistics are then used to quantify the magnitude of deviation from the normative model using a distribution of extreme Z-scores. Then, deviations falling within a specified range (e.g. in the top or bottom 1%) of the extreme Z-score distribution are to be considered abnormal relative to the normative model.

Here abnormality was defined by a critical Mahalanobis distance value computed from Wilk’s criterion based on conservative pre-defined confidence level. Given the nature of the injuries this approach is appropriate since expected magnitude of changes is large in the patients that were studied. However, there may be other conditions, such as mild TBI or autism, where individual variations may be more subtle. Thus, depending on the application, a reduced penalization for multiple comparisons may be desired so that the threshold for abnormality could be made less conservative. In that case, an approach such as that in Marquand et al. (2016) may be adopted for empirical determination of a subjectlevel abnormality index if there are enough participants in the reference group. For example, a normal range for the Mahalanobis distance could be defined as the 90–95% confidence interval based on the distribution of the control data.

Limitations and future work

This study has a number of limitations that should be considered when interpreting the results. First, it is well known that the interpretation of DTI is ambiguous in crossing fiber regions. The MaD-Tract method itself is agnostic to crossing fibers since it focuses on developing a distribution of eigenvalues at each location along a tract and, if significantly different from the normative distribution, distributions are flagged as abnormal. Nonetheless, the technique could certainly be less sensitive to changes in areas of crossing fibers or complex fiber organization and that warrants further exploration. In future work, the nique could be extended to include higher order models like NODDI, which directly address fibers dispersion within a voxel. Second, measurement variation across scanners and sites cannot be entirely ruled out as a contributing factor in this study. Normative control data and/or human traveling phantoms at each site would have been ideal, but unfortunately were not feasible. However, measurements with a standardized PVP diffusion phantom across several of the sites revealed remarkably consistent diffusion estimates across those sites, particularly with respect to Site I where the control participants were scanned. Model-based correction of the estimated eigenvalues using the phantom data did not alter the abnormal findings of those cases. Third, while the selection of a critical MaD value for abnormality accounts for the reference sample size, the method’s performance should improve with better estimates of the covariance matrix from increased size of the normative control group. Further, in this application, age effects in the DTI eigenvalues were not accounted for as these were observed to be small effects. Nonetheless, age dependence does impact the characteristics of the normative group distribution by adding to the variance and modeling out those effects may improve the sensitivity of the method. Additionally, this method is currently limited to white matter regions, but, as is the case with other neurological conditions, traumatic lesions in TBI include both white and gray matter (Avants et al., 2008; Kurča et al., 2006). Finally, Mad-Tract results are not discussed in the context of clinical outcomes of severe pediatric TBI. While the purpose of this work was to stablish the framework, the relations that may exist between this multivariate score and the development of longer-term complications needs to be investigated. The analyzed tracts were selected to ultimately relate brain injuries to clinical impairment and global cognitive function (corpus callosum), executive function (cingulum bundles), and mood and affective function (uncinate fasciculus). This is the subject of ongoing work.

Supplement

For the subjects scanned at sites different from the control group the figures below showcase the skeletonized tracts with detected abnormalities embedded in the subjects T1w volumes.

Patient C

Patient D

Patient M

Patient R

Patient S

Patient V