Dimensionality reduction of diffusion MRI measures for improved tractometry of the human brain.

Various diffusion MRI (dMRI) measures have been proposed for characterising tissue microstructure over the last 15 years. Despite the growing number of experiments using different dMRI measures in assessments of white matter, there has been limited work on: 1) examining their covariance along specific pathways; and on 2) combining these different measures to study tissue microstructure. Indeed, it quickly becomes intractable for existing analysis pipelines to process multiple measurements at each voxel and at each vertex forming a streamline, highlighting the need for new ways to visualise or analyse such high-dimensional data. In a sample of 36 typically developing children aged 8-18 years, we profiled various commonly used dMRI measures across 22 brain pathways. Using a data-reduction approach, we identified two biologically-interpretable components that capture 80% of the variance in these dMRI measures. The first derived component captures properties related to hindrance and restriction in tissue microstructure, while the second component reflects characteristics related to tissue complexity and orientational dispersion. We then demonstrate that the components generated by this approach preserve the biological relevance of the original measurements by showing age-related effects across developmentally sensitive pathways. In summary, our findings demonstrate that dMRI analyses can benefit from dimensionality reduction techniques, to help disentangling the neurobiological underpinnings of white matter organisation.

Introduction

The human brain is composed of multiple white matter fibres connecting gray matter areas dedicated to processes such as memory, cognition, language, or consciousness. Diffusion MRI (dMRI) (Basser et al., 1994, 2000; Basser and Jones, 2002; LeBihan et al., 2001) has become the preferred tool to probe the brain's tissue microstructure non-invasively. Measures derived from diffusion tensor imaging (DTI) (Basser et al., 1994) can be obtained at each imaging voxel, including fractional anisotropy (FA) which reflects the degree of diffusion anisotropy (Pierpaoli and Basser, 1996), and mean diffusivity (MD), an indicator of the overall magnitude of diffusion. Based on local estimates of underlying trajectories at every voxel, dMRI is also capable of virtually reconstructing the structural architecture of the brain white matter pathways using tractography (Conturo et al., 1999; Mori and Van Zijl, 2002). The conventional approach to merge the quantitative nature of diffusion measures with the qualitative nature of tractography is to collapse voxel-based measures into a single scalar value per bundle (e.g., by averaging values over all vertices of a streamline; Jones et al. (2006); Kanaan et al. (2006); Jones et al. (2005a)). Individual differences in such summary diffusion-related measures can then be correlated, for example, with individual differences in cognition or behaviour. However, despite its well documented sensitivity, DTI has its limitations (Tournier et al., 2011; Jeurissen et al., 2013). For example, FA and MD lack specificity to the various physical properties of white matter, such as crossing fibres (Jeurissen et al., 2013), axon density and myelination (Beaulieu, 2002; Jones et al., 2013). Moreover, the average profile of those measures may vary along a given pathway depending on the underlying fibre architecture (Vos et al., 2012; Yeatman et al., 2012). Furthermore, only a subset of DTI measures are known to be orthogonal with each other (e.g., FA, MD or tensor norm; Ennis and Kindlmann (2006); Kindlmann et al. (2007); De Santis et al. (2014)).

Recent advances in diffusion hardware, acquisition and modelling (Sotiropoulos et al., 2013; Jones et al., 2018; Assaf and Basser, 2005; Tournier et al., 2012; Jeurissen et al., 2014) have been introduced to overcome the limitations of DTI, giving access to previously inaccessible measures. High angular resolution diffusion imaging (HARDI; Tuch et al. (2002)) was originally developed to not only provide new anisotropy measures (Tournier et al., 2011) but also to solve the so-called crossing fibre problem, making tractography more robust (Descoteaux, 2015). Multi-shell acquisitions (Wedeen et al., 2005) have also facilitated new ways to link relevant tissue properties to the signal such as CHARMED (Assaf and Basser, 2005), AxCaliber (Assaf et al., 2008), ActiveAx (Alexander et al., 2010), multi-tensor models (Scherrer et al., 2016) and NODDI (Zhang et al., 2012) among others (for review, see Alexander et al. (2017)). In general, such models aim to extract parameters from intra- and extracellular compartments, and to estimate parameters such as axon diameter distributions and other high-order information.

Multi-shell acquisitions have also shown to improve the angular resolution of orientation distribution functions (ODFs) (Descoteaux et al., 2011; Jeurissen et al., 2014; Chamberland et al., 2018). In conjunction, new frameworks such as fixel-based analysis (Raffelt et al., 2012) have been proposed to map fibre-specific measures by looking at the apparent fibre density (AFD), a measure proportional to the underlying fibre density, as opposed to having voxel-specific scalar maps. The combination of frameworks such as along-tract profiling (Jones et al., 2005b; Corouge et al., 2006; Yeatman et al., 2012; De Santis et al., 2014; Colby et al., 2012; Cousineau et al., 2017) and tractometry (e.g., combining multiple measures (Bells et al., 2011)) allows for a comprehensive assessment of white matter microstructure. Both frameworks have the advantage of providing higher sensitivity to microstructural features of fibre pathways by mapping a set of MR-derived measures over white matter bundles. Recently, along-tract profiling has been successfully applied to study normal brain development (Geeraert et al., 2018) and to characterise areas of the brain with abnormal properties in various brain conditions (Dayan et al., 2016; Cousineau et al., 2017; Groeschel et al., 2014).

However, one problem arises with having access to multiple new measurements at each voxel and at each vertex forming a streamline: it quickly becomes intractable for existing analysis pipelines to process such high-dimensional data (a problem often referred to as the curse of dimensionality; Bellman (1961)), highlighting the need for new ways to visualise or analyse such data. Moreover, dMRI measures may share overlapping information which can cause redundancies (in the sense of correlation) in data analysis and ultimately decrease statistical power if strictly correcting for Type I errors (Penke et al., 2010; Metzler-Baddeley et al., 2017; Bourbon-Teles et al., 2017). A solution to this problem resides in dimensionality reduction, an established technique that has been successfully applied in the past by the neuroimaging community (for review, see Mwangi et al. (2014)). Despite the growing number of experiments using different microstructural measures in assessments of white matter, there has been limited work on combining these different measures and on examining their covariance along specific pathways.

In this work, we explore the covariance of commonly-derived dMRI measures (De Santis et al., 2014). We propose a data reduction framework that takes advantage of those redundancies and aims to provide a better insight into patterns of associations between DTI and HARDI measures. Specifically, we identified common components that explain the maximal variance in measures profiled along multiple fibre bundles. We demonstrate the utility of our framework by showing enhanced sensitivity to the detection of age-related differences in tissue microstructure across developmentally sensitive pathways compared with the individual dMRI measures. Finally, we provide recommendations for future studies with limited capabilities in terms of data acquisition and processing.

Methods

Participants

This study reports on a sample of typically developing children aged 8–18 years (mean = 12.2 ± 2.8) participating in the Cardiff University Brain Research Imaging Centre (CUBRIC, School of Psychology) Kids study. The study was performed with ethics approval from the internal ethics review board and informed consent was provided from the primary caregiver of children enrolled in the study. Exclusion criteria included previous history of a neurological condition or epilepsy.

Data acquisition

Data from thirty-six (n = 36, 13 male) children were acquired using a multi-shell HARDI protocol on a Siemens 3T Connectom system with maximum gradient amplitude = 300 mT/m. The acquisition protocol consisted of 14 b0 images, 30 diffusion directions at b = 500, 1200 s/mm2 and 60 diffusion directions at b = 2400, 4000, 6000 s/mm2 with 2 × 2 × 2 mm3 voxels (TE/TR: 59/3000 ms, δ/Δ: 7.0/23.3 ms).

Data pre-processing

Data quality assurance was performed on the raw diffusion volumes using slicewise outlier detection (SOLID; Sairanen et al. (2018)). Each dataset was then denoised in MRtrix (Veraart et al., 2016) and corrected for signal drift (Vos et al., 2017), subject motion (Andersson and Sotiropoulos, 2016), field distortion (Andersson et al., 2003), gradient non-linearities (Glasser et al., 2013; Suryanarayana et al., 2018) and Gibbs ringing artefacts (Kellner et al., 2016).

Local representation

Multi-shell multi-tissue constrained spherical deconvolution (MSMT-CSD; Jeurissen et al. (2014)) was applied to the pre-processed images to obtain voxel-wise estimates of fibre ODFs (fODFs; Tournier et al. (2004, 2007); Seunarine and Alexander (2009); Descoteaux et al. (2009)) with maximal spherical harmonics order  = 8. The fODFs were generated using a set of 3-tissue group-averaged response functions (Dhollander et al., 2016) followed by joint bias field and image intensity normalisation in MRtrix (Tournier et al., 2012), enabling the direct comparison of fODF amplitudes across subjects (Raffelt et al., 2012). Diffusion tensors were also generated using linearly weighted least squares estimation (for b < 1200 s/mm2 data) providing the following quantitative scalar measures: FA, axial diffusivity (AD), radial diffusivity (RD), MD, geodesic anisotropy (GA; Fletcher et al. (2004)) and tensor mode representing the shape of the tensor (Kindlmann et al., 2007). In addition, HARDI measures were extracted from the fODFs of each subject. Those measures include fibre-specific AFD (Raffelt et al., 2012) for the bundles described in the next section, AFD (spherical harmonics l = 0) and the Number of Fibre Orientations (NuFO) based on the number of local fODF peaks (Dell’Acqua et al., 2013). Finally, restricted signal fraction maps (FR, adapted from CHARMED to remove potential isotropic partial volume contamination; Assaf and Basser (2005)) were also computed using the fODFs peaks to initialise and regularise model-fitting. To summarise, ten dMRI measures related to tissue microstructure (m = 10) were generated for each subject.

Tractography and tractometry

Whole-brain streamline tractography was performed using FiberNavigator (Chamberland et al., 2014) using 8 seeds/voxel evenly distributed across the whole brain (approximating 1.8 M seeds), a minimum fODF amplitude of 0.1, a 1 mm step size (i.e. 0.5voxel size), a 45∘ maximum curvature angle and streamlines whose lengths were outside a range of 20 mm–300 mm were discarded. Twenty-two bundles of interest (t = 22) were then interactively dissected in the native space of each subject using a combination of include and exclude regions of interests (ROIs). Anatomical definitions and ROIs used to delineate each pathway are listed in the Supplementary Materials. The virtual dissection plan included:

Commissural bundles: anterior commissure (AC), body of the corpus callossum (CC), forceps minor (Genu), forceps major (Splenium).

Association bundles (bilateral): arcuate fasciculus (AF), cingulum (Cg), inferior fronto-occipital fasciculus (iFOF), inferior longitudinal fasciculus (ILF), optic radiations (OR), superior longitudinal fasciculus (SLF), uncinate fasciculus (UF).

Projection bundles (bilateral): corticospinal tract (CST), frontal aslant tract (FAT).

At this stage, we examined the covariance of the averaged diffusion measures for all bundles using Pearson's correlation (r). Next, along-tract profiling was performed for each bundle using the Python toolbox developed by Cousineau et al. (2017), combined with DIPY functions (Garyfallidis et al., 2014). Bundles were first pruned to remove outliers as in (Cousineau et al., 2017; Garyfallidis et al., 2018). If necessary, the order in which the vertex-wise measures were stored was reversed, to ensure consistency across the subjects profiles (i.e., from left-to-right for commissural bundles, from inferior-to-superior for projection bundles and from posterior-to-anterior for association bundles). A representative core streamline was generated for each bundle (i.e., mean streamline of the pathway) and was subsequently resampled to s = 20 equidistant segments. Then, every vertex of every streamline forming the pathway was assigned to its closest segment along the core. The measure values of each vertex were then projected and averaged along each segment of the pathway, weighted by their geodesic distance from the core (Cousineau et al., 2017). An along-tract profile was finally generated for every combination of measure and pathway.

Dimensionality reduction

Each dataset comprised m = 10 dMRI-derived measures mapped along 440 white matter regions (t = 22 bundles × s = 20 segments). To explore the possible redundancy (in the context of data reduction) and complementarity of each measure, a principal component analysis (PCA) was performed on the concatenated set of profiles across subjects and bundles (Table 1) using the tidy data standard (Wickham et al., 2014). Performing PCA over tract segments rather than over all voxels alleviates the need to register each diffusion measure to a common space. PCA reduces data dimensionality by extracting principal components that reflect relevant features in the data (Jolliffe, 2002; Abdi and Williams, 2010). The benefit is that a significant proportion of the variance in the data can be explained by a reduced number of orthogonal components, compared to the total number of raw input variables. PCA was performed by singular value decomposition of the z-transformed tract profiles via the prcomp function in R (R Core Team, 2018). Here, the goal was to end up with the minimum number of components that summarise the maximum amount of information contained in the original set of diffusion measures. However, in order to avoid instability around the component loadings that comprise the principal components (Garg and Tai, 2013), measures showing significant covariance were discarded based on their correlation scores (|r| > 0.8) and the PCA was re-computed. Finally, the minimal number of principal components that accounted for the most variability was selected based on: 1) their interpretability (Metzler-Baddeley et al., 2017); and 2) the inspection of scree plots (Cattell, 1966) to select ranked components with an eigenvalue > 1.

Table 1

Data structure input for PCA. Individual subjects (n = 36), bundles (t = 22) and segments (s = 20) are concatenated to form observations while variables represent the measures (m = 10) derived from dMRI.

Subject	Bundle	Section	FA	AD	…	FR
S1	Bundle1	Section1	FA111	AD111	…	FR111
S2	Bundle1	Section1	FA211	AD211	…	FR211
⋮	⋮	⋮	⋮	⋮	…	⋮
S1	Bundle1	Section2	FA112	AD112	…	FR112
⋮	⋮	⋮	⋮	⋮	⋱	⋮
Sn	Bundleb	Sections	FAnbs	ADnbs	…	FRnbs

Statistical analysis

PCA results were tested for sampling adequacy using a Kaiser-Meyer-Olkin (KMO; Dziuban and Shirkey (1974)) test followed by Bartlett's test of sphericity to test whether the covariance matrix is significantly different from identity. We then ran an exploratory linear regression analysis to see whether profiles extracted from the PCA can provide increased sensitivity in the detection of age-related differences in tissue microstructure (as opposed to using the full set of m = 10 measures). It is important to recall that PCA results are always orthogonal, and therefore are statistically independent of one another. To address the multiple comparisons problem, a strict Bonferroni correction was applied to all linear models whereby statistical significance was defined as: p < 0.05/(m measures × t bundles × s segments) resulting in p < 1.14e-5 for the ten raw measures, and p < 5.68e-5 for the first two principal components. All statistical analyses were carried out using R v3.5.1 (R Core Team, 2018) and RStudio v1.1.456 (RStudio Team, 2016).

Results

Measures covariance and profiling

The entire set of bundles and measures was successfully reconstructed in all subjects. Fig. 1 shows the relationship between the various input measures averaged on different white matter pathways using a cross-correlation matrix representation. Matrices are re-organised using hierarchical clustering (Murtagh and Legendre, 2014), placing higher correlations closer to the diagonal in order to regroup measures that have similar correlations together. In general, the measures form two or three clusters across all bundles. Across the whole set of bundles (Fig. 1, middle), we observe a first cluster of positive correlations (r > 0.5) between AD, FA, GA, Mode, AFD, AFD and FR measures. A second cluster of positive correlations is formed of MD, RD and NuFO. The group-averaged diffusion measures of each bundle are reported in Suppl. Table 1.

Fig. 1

Correlation matrices of the ten diffusion measures, group-averaged for each extracted bundles. The middle image represents the average of all white matter bundles. Matrices are re-organised using hierarchical clustering, grouping measures that have similar correlations together. A first cluster of positive correlations (r > 0.5) is observed between most of the bundles for measures like AD, FA, GA, AFD, AFD, Mode and FR. A second set of positively correlated measures (NuFO, MD, RD) forms the second cluster. Note that for bilateral pathways, the left and right values were combined prior performing the correlation.

However, important details about the spatial heterogeneity of the various input measures of interest appear when profiled along pathways. For example, FA, AFD and FR values all get progressively smaller along the CST as they approach the cortex (Fig. 2, Fig. 3). Furthermore, the high number of fibre crossings near the centrum semiovale is reflected by a high NuFO index and is also marked by a low FA (Fig. 2). As might be expected, HARDI-derived measures such as FR and AFD seem to be less affected by the intra-voxel orientational heterogeneity of crossing regions than the tensor-based measures like FA, AD and RD. The correlation matrix in Fig. 2, Fig. 3 also highlights the similarity between the various microstructural profiles, indicating potential overlap in the amount of information conveyed by the dMRI measures.

Fig. 2

Overview of the ten input measures overlaid on the CST of a representative subject. Whole-brain tractograms (top-left) were manually dissected into t = 22 bundles (bottom-left) and measures were subsequently mapped along each pathway, providing information about their spatial heterogeneity.

Fig. 3

Group-average profiling of the ten input diffusion measures along the left CST for s = 20 segments, spanning from the brainstem (s = 1) to the cortex (s = 20). Heterogeneity in the profiles along the tract highlights the need for a vertex-wise assessment of the measures. Similarity between profiles also shows shared covariance between the measures, indicated by the two clusters (1 and 2) on the correlation matrix (sorting: hierarchical clustering). Shaded tract-profile area: ±1 standard deviation.

Principal component analysis

The loading vector plots in Fig. 2, Fig. 3 show association patterns between the various input measures. The left panel shows PCA results performed on the entire set of measures. If two vectors subtend a small angle to each other, the two variables they represent are strongly correlated. When such vectors were found to be close, the one showing higher correlations with any other measures was removed (see Section 2.6). In line with the aforementioned results, shared covariance is observed between AD and tensor mode (r = 0.8), as well as between FA and GA (r = 0.95). After pruning up measures for multicollinearity (Fig. 2, Fig. 3, right panel), PCA results show that 80% of the variability in the data is accounted by the first two principal components (KMO: 0.64, sphericity: p < 2.2e-16). As shown in Fig. 2, Fig. 3, the PC that explains the largest proportion of the variance (PC1, 48%, λ = 3.4) is composed of hindrance-sensitive measures (Dell’Acqua et al., 2013) with AFD, FR and AFD contributing positively (24%, 21% and 16%, respectively) and AD contributing negatively (25%). The second PC (PC2) represents 32% of the variance in the data (λ = 2.2) and is mostly driven by orientational dispersion and complexity-sensitive measurements, with its largest positive contribution from NuFO (34%), and negative contributions from AD (26%) and MD (25%) Fig. 2, Fig. 3, PC2).

PCA results before (left) and after (right) multicollinearity analysis. To improve stability around PC1, AFD was kept over GA and FA due to its fibre specificity properties. Tensor mode was also discarded based on its collinearity with AD. On the right PCA, one can observe separation between the various measures, generating a hindrance-related component (PC1) that loads on AFD, AFD, RD and FR and a complexity-related component (PC2) that loads on NuFO, AD and MD. Here, the squared cosine notation (cos2) shows the importance of a measure for a given PC. A high cos2 indicates a good representation of the measure for a given principal component. Figure generated using the FactoExtra package (Kassambara and Mundt, 2017) in R.

Visual overview of PC1 and PC2 with the contribution of each measure to that component. The first PC captures most of the hindrance- and restriction-related measures (AFD, RD, FR, AFD). The second PC mostly represents tissue complexity and orientational dispersion properties associated with NuFO, AD and MD. The first two components accounted for 80% of the variance in the diffusion measures.

Detecting differences in tissue microstructure

Bundle averages

To asses the relevance of the two principal components, we first report developmental changes in white matter tissue microstructure using bundle-averaged measures (m = 10, s = 1) and PCA components (m = 2, s = 1). Significance thresholds were Bonferroni-corrected to account for multiple comparisons (p < 2.27e-4 for the ten raw measures, and p < 1.14e-3 for the two principal components). Fig. 2, Fig. 3 shows a significant increase in PC1 as a function of age for the left iFOF and CST, whereas no correlation with age was observed between individual hindrance-related measures. Significant positive correlations were found between PC1 (restriction-related component) and age in the following association pathways: AF (left: R2: 0.34, p = 1.06e-4), iFOF (left: R2: 0.31, p = 2.51e-4), FAT (left: R2: 0.43, p = 9.06e-6, right: R2: 0.43, p = 7.77e-6), UF (right: R2: 0.26, p = 9.76e-4) and motor pathways: CST (left: R2: 0.40, p = 1.92e-5, right: R2: 0.40, p = 2.0e-5), CC (R2: 0.29, p = 4.51e-4). One significant positive correlation between PC2 (complexity-related component) and age was found in the SLF (right: R2: 0.27, p = 6.82e-4).

Developmental changes in iFOF and CST bundles. PC1 show significant positive correlation with age, whereas no correlation was observed between the individual hindrance-sensitive measures.

No significant age relationships were found between any of the bundles and FA, GA, Mode, AFD, AFD, FR and NuFO. Significant negative correlations with age were found in RD for the AF (left: R2: 0.37, p = 5.40e-5, right: R2: 0.34, p = 1.02e-4), Cg (left: R2: 0.33, p = 1.45e-4, right: R2: 0.47, p = 2.22e-6), CST (left: R2: 0.35, p = 7.36e-5, right: R2: 0.33, p = 1.28e-4), FAT (left: R2: 0.37, p = 4.65e-5, right: R2: 0.28, p = 4.67e-4). Significant negative correlations with age were found in MD for the FAT (left: R2: 0.33, p = 1.22e-4, right: R2: 0.31, p = 2.18e-4), AF (left: R2: 0.39, p = 2.48e-5, right: R2: 0.36, p = 6.28e-5), CST (left: R2: 0.40, p = 2.06e-5, right: R2: 0.42, p = 1.12e-5), SLF (right: R2: 0.38, p = 3.21e-5) and Cg (right: R2: 0.37, p = 4.06e-5). One significant negative correlation was found between AD and age in the SLF (right: R2: 0.41, p = 1.80e-5).

Along-tract profiling

Here, we report on developmental changes in tissue microstructure seen with along-tract profiling (see Fig. 7). Table 2 reports the measures and tract segment mapped along different pathways where significant correlation with age was observed. Significant positive correlations were found between PC1 (restriction-related component) and age near the motor cortex area for the CST20 (right: R2: 0.37, p = 4.91e-5) and CC3 (R2: 0.37, p = 4.62e-5). Significant age-related positive correlations with PC2 were observed for motor-related pathways: CC2 (R2: 0.37, p = 5.63e-5), CST20 (right: R2: 0.49, p = 1.41e-6), CST19 (left: R2: 0.43, p = 8.59e-6) and association pathways: FAT19 (right: R2: 0.38, p = 4.08e-5), iFOF19 (right: R2: 0.39, p = 2.89e-5), ILF1 (right: R2: 0.35, p = 9.26e-5), UF3,4 (right: R2: 0.43 & 0.39, p = 7.71e-6 & 3.0e-5) and SLF1,2 (R2: 0.46 & 0.42, p = 3.41e-6 & 1.14e-5). Suppl. Fig. 2 shows that the PCA profiles preserve the spatial heterogeneity of the input diffusion measures along the CST.

Fig. 7

Age relationships captured by PC1 and PC2 over the left CST. Highlighted section of an axial slice overlaid with fODFs reconstruction of a representative participant (top left) shows the contoured area (black line) where streamlines terminated to form segment 20. At the group level, significant positive correlations with age were found with PC1 and PC2 (top right). Significant positive correlations were also found for HARDI measures AFD and NuFO (middle right). No significant correlations were observed for any of the DTI measures (bottom right). Profile plots indicate where significant differences in tissue microstructure were located along the CST (-log(p) scale).

Table 2

Segments of white matter bundles where significant correlation between diffusion measures and age was observed. Subscript ordering for along-tract positions: left (s = 1) to right (s = 20) for commissural bundles, inferior (s = 1) to superior (s = 20) for projections bundles and posterior (s = 1) to anterior (s = 20) for associations bundles. Positive and negative correlations are indicated by increasing (↗) and decreasing (↗) arrows, respectively. Significance thresholds for the measures and components were set as p < 1.13e-5 and p < 5.68e-5, respectively (adjusted R2 > 0.3).

Individual diffusion measures					PCA
RD	MD	AFDtot	NuFO	FR	PC1	PC2
↘ r-FAT8	↘ r-AF16	↗ l-CST1,20	↗ CC1	↗ r-FAT6	↗CC3	↗CC2
	↘ l-Cg7	↗ r-CST20	↗r-CST20		↗r-CST20	↗r-CST20
						↗ l-CST19
						↗ r-FAT19
						↗ r-iFOF19
						↗ r-ILF1
						↗ r-UF3,4
						↗ r-SLF1,2

No significant age relationships were found in any bundles for FA, GA, Mode, AD and AFD. Significant negative correlations were observed for RD in the FAT8 (right: R2: 0.42, p = 1.08e-5) and for MD in the AF16 (right: R2: 0.44, p = 5.62e-6) and Cg7 (left: R2: 0.43, p = 6.5e-6). Significant positive correlations were also found for AFD in the CST20 (left: R2: 0.50, p = 9.02e-7, right: R2: 0.53, p = 2.73e-7) and CST1 (left: R2: 0.46, p = 2.88e-6). For NuFO, significant age-related differences in tissue complexity were observed in the CST20 (right: R2: 0.54, p = 1.56e-7) and CC3 (R2: 0.43, p = 7.09e-6). Finally, one significant positive correlation in FR was found for the FAT6 (right: R2: 0.57, p = 4.9e-8).

Discussion

Extraction of interpretable components

The aim of this study was to systematically examine any potential covariance between various diffusion measures mapped along white matter fibre bundles extracted from a cohort of typically-developing children and adolescents. We first examined the covariance of the measures averaged over different bundles, which revealed two clusters of inter-dependent measures. The first cluster revealed that measures sensitive to restricted diffusion shared high correlations with each other. Similarly, measures which are known to be sensitive to local complexity or orientational-dispersion co-varied. When profiled along pathways, measures showed heterogeneity across the trajectory of diverse pathways, but their interdependency remained marked by a two-cluster formation. This provided motivation for the next step of our analysis, where we performed a PCA to collapse the inter-dependent measures into the principal modes of variance. This was done by profiling multiple brain fibre systems based on their dMRI features, and then deriving a set of principal components that best represent those individual measures. We then showed the sensitivity of these new components to the detection of differences in tissue microstructure of white matter pathways by exploring their relationship with the age of participants.

A common problem with PCA is that the interpretation of the resulting components can be challenging. Here, the principal components loaded onto variables that shared similarities in their sensitivity to different tissue properties, making the interpretation of the resulting components more meaningful. Measures accounting for the largest percentage of variance in the data (forming PC1) are those known to be most sensitive to hindrance or restriction in the signal, including RD, AFD and FR. In contrast, PC2 features measures that could reflect complexity or orientational-dispersion in the signal, such as NuFO, AD and MD.

The raw tract-averaged diffusion measures showed significant negative correlations in MD and RD with age across a range of developmentally sensitive tracts, which is in line with previous studies that also report a decrease in MD and RD with age, whereas FA shows slower increases with age in late childhood (Lebel et al., 2008). Additionally, we observed significant age relationships for PC1 and PC2 in developmental-sensitive tracts involved in language (i.e., AF, SLF, FAT, iFOF, ILF) and motor functions (i.e., CST and CC), which is also in line with previous reports in the literature (Genc et al., 2017, 2018, 2018; Lebel et al., 2017; Geeraert et al., 2018). Examining the tract-averaged values for PC2, only one significant correlation with age was found for the right SLF. However, additional significant correlation between age and PC2 were observed when performing along-tract profiling. A potential explanation for this findings resides in the nature of PC2, with most of its contributions coming from AD, MD and NuFO. Indeed, since PC2 reflects the local complexity at each voxel, measures like NuFO will vary depending on the underlying structural architecture and therefore taking the average value across the bundle may lead to a summary statistic that is hard to interpret. In contrast, the values of PC1 (or AFD) remain relatively constant over bundles and thus are less impacted by calculating the bundle average. This may also explain why the first principal component derived from bundle-averages was found to show significant correlations with age in white matter bundles, whereas the original DTI measures did not. We also note that sex differences were not accounted for due to our relatively small sample size but may play a role in some of the differences we observed across pathways (Seunarine et al., 2016). Moreover, studies of larger size would be better powered to test whether the resulting PCs are linked with other factors over and above the participant's age, such as sex, IQ and cognition. Yet, our results highlight the sensitivity of PC1 and PC2 as a composite measure by 1) showing significant correlation with age in regions where other measures did not, and 2) reflecting effects captured by the other measures.

Restricted (or hindered) diffusion is primarily caused by dense packing of axons and their cell membranes (Beaulieu, 2002). Other tissue properties such as myelination and local complexity can also affect the degree of hindrance or restrictance measured at each voxel (Vos et al., 2012). In the current study, an increase in PC1 may indicate higher coherence of the underlying white matter bundles for our older subjects, in comparison with the younger ones. Previous studies have demonstrated that dMRI measures can be sensitive to age-related differences, and those are often associated with an increased microstructural organisation (for review, see Lebel et al. (2017)). Given the well-established role of the CST in supporting motor performance, our finding of increased hindrance with age in typically developing children is in line with previous research that showed that brain maturation varies across different pathways, with commissural and projection tracts reaching maturation by early adolescence while association pathways develop over a longer time period (Geeraert et al., 2018). Interestingly, PC2 also captured an increase in orientational dispersion for that same region (which was either marked by an increase in NuFO or decrease in MD, Fig. 2, Fig. 3). The fact that those changes appear near the cortex, a region usually contaminated by partial volume effects, highlights the role of MSMT-CSD in achieving adequate fODFs representation at the boundary between gray matter and white matter (Jeurissen et al., 2014).

In light of our results, we stress that the proposed framework was applied to study neuro-developmental changes related to age only, and therefore, results should be interpreted with care in any other context. The framework should be used as a general approach to data reduction of MRI-derived measures. However, this does not suggest that all studies will benefit from this form of data reduction since our results may not be generalizable to other populations. Thus, we recommend that our proposed framework is applied to each new data set to discover the most appropriate PCA loadings. Indeed, it is entirely plausible that in a disease population or data collected with different acquisition parameters or hardware, the optimal factor loadings would be different than the ones found here.

Choice of diffusion measures

A growing interest in utilising advanced dMRI measures to study the human brain motivated us to investigate the shared relationship between DTI and HARDI measures (De Santis et al., 2014). Ultimately, the signals captured from a white matter bundle are coming from the same substrate i.e., from the same axons, myelin and other cellular inclusions, so some degree of correlation is likely to be observed. Redundancy between the different measures (in the sense of correlation) does not however imply that the some of the measures are not useful. In the event that two measures are strongly correlated, any deviation from perfect correlation may reflect that each is capturing subtly different information that may indeed be crucial (for example in understanding disease processes) and thus, regarding one redundant measure to be discarded in favour of the other may not be sound decision. Here, we adopt the term redundancy to refer explicitly to data reduction only.

Being a relatively fast-developing field, dMRI offers a multitude of mathematical models to represent the underlying tissue microstructure (Alexander et al., 2017). Here, we focused on DTI and HARDI measures (Descoteaux, 2015). dMRI measures are generally sensitive to differences in tissue microstructure that can potentially be linked to fibre properties such as myelination and axon density (Scholz et al., 2014). Despite the fact that the specific interpretation of these measures remains controversial (Jones et al., 2013), DTI and HARDI measures are routinely used by neuroscientists and clinicians to gain insights into white matter properties. The findings reported here are in line with existing evidence suggesting that HARDI measures may be more specific than DTI for the detection of differences in tissue microstructure (Tournier et al., 2011; Jeurissen et al., 2013; Cousineau et al., 2017). Our results suggest that combining the sensitivity of DTI and the specificity of HARDI has the potential for compromise between the two techniques. Other macrostructural measures (e.g., bundle volume or mean length; Lebel et al. (2012, 2008, 2017); Geeraert et al. (2018); Girard et al. (2014)) have been used to study brain development and may also provide complementary features that could be added in the proposed framework. Other measures such as rotationally invariant spherical harmonic features (RISH; Mirzaalian et al. (2015); Caruyer and Verma (2015); Zucchelli et al. (2018)) could also be introduced in the current framework, with the main advantage of representing more directly the diffusion signals rather than relying on various microstructural models. Ultimately, the key challenge resides in knowing what measure (or combination of; De Santis et al. (2014)) provides the best value in terms of scanning and processing time. To help with the planning of future studies and based on our observations, we present some recommendations for data analysis.

DTI and HARDI: DTI measures can nowadays be easily be derived from a conventional 30 directions protocol acquired at b = 1000 s/mm2 in approximatively 5 min. Here, all DTI measures were derived using only the b = 1200 s/mm2. HARDI measures such as AFD (Raffelt et al., 2012; Dell’Acqua et al., 2013) and NuFO (Dell’Acqua et al., 2013) are also readily derived using a standard 3T MRI scanner. Indeed, CSD can usually be performed on data acquired with a minimum of 30–45 directions at b = 1000 s/mm2 (Dell’Acqua and Tournier, 2018; Alexander and Barker, 2005; Tournier et al., 2013) or even 21 directions if the quality is acceptable (Chenot et al., 2018a). Moreover, going beyond single b-value acquisitions will provide a better estimation of partial volume effects and better characterisation of various tissue types that will subsequently improve HARDI reconstructions in those areas (Jeurissen et al., 2014; Chamberland et al., 2018).

Tract-averaging and along-tract profiling: In the context of along-tract profiling, age-related effects might be more pronounced when performing group-wise comparisons (e.g., young vs old, patients vs controls; Yeatman et al. (2014)) rather than directly looking at a single cross-sectional change in tissue microstructure. Indeed, theses changes might be too subtle to detect, especially considering that the age range of our participants falls on the inflection point of the developmental curve (Lebel and Beaulieu, 2011). Moreover, one may consider first looking at the profile along each tract of interest and ask the following: are there any benefits in sub-segmenting the profile into finer portions? Admittedly, if the measure of interest remains stable along the pathways, a conventional tract average is probably better suited than looking at a constant profile at multiple points. Depending on the research hypothesis, the use of a more permissive approach such as false discovery rate (FDR) correction could be considered to assess differences along multiple adjacent bundle segments. Another possible approach to analyse the multi-dimensionality of the data resides in functional data analysis (FDA; Ramsay (2005); Ferraty and Vieu (2006); Wang et al. (2016)), a statistical method that operates on continuous or discrete functions (e.g., tract profiles; Goodlett et al. (2009)) and that takes into account the spatial interdependency of each segment. In addition, the cluster patterns shown in Fig. 1 may suggest that performing a PCA on each bundle separately might result in different decompositions. We presume this approach would be more appropriate in the context of a bundle-specific analysis (or hypothesis-driven approach). However, by doing so, comparisons between bundles become impossible, as there is no guarantee that e.g., the first PC of a given bundle could be matched with the first PC of another bundle.

Lastly, we acknowledge that the data acquisition and processing employed in this study were performed on unique hardware; however, this should not discourage future users to adopt the current framework for their own data analysis. Again, one should not infer that we are recommending that the factor loadings presented here should be used off-the-shelf; rather, we are recommending that the general framework presented here is applied to the specific application, acquisition, and hardware so that the data reduction is tailored. Yet, advances in data harmonisation (Tax et al., 2019) show great promise in bridging the gap between state-of-the-art and standard data acquisition schemes.

Future perspectives

Over the last few years, a wide range of supervised and unsupervised learning applications based on feature extraction has emerged, ranging from individual classifiers for specific brain disorders (Wang et al., 2010; Chu et al., 2012) to data predictors of brain function (Chen et al., 2009; Franke et al., 2012; Casanova et al., 2011; Kucukboyaci et al., 2014). In the field of functional MRI, independent component analysis (ICA) is a successful example of unsupervised dimensionality reduction that allows the extraction of temporally segregated resting-state networks (Beckmann et al., 2009). In all cases, data reduction approaches facilitate a stream to analyse and interpret the increasingly large multi-dimensional data generated by new methodological models. Admittedly, despite PCA being one of the most commonly-used tool for data reduction, it can also over-fit data (Kramer, 1991) and potentially require multiple post-processing regression steps to explore the link between the resulting components and the observations. Similar to PCA, a canonical correlation analysis (Hotelling, 1936) may help in finding the link between the correlated measures and observations by extracting their joint information. Other non-linear dimensionality reduction techniques based on manifold learning such as isometric feature mapping (Tenenbaum, 1998) or locally linear embedding (Roweis and Saul, 2000) may also better disentangle the measurement space. However, one has to be cautious when applying advanced models of dimensionality reduction to medical imaging as it is often a trade-off between model accuracy and interpretability. Indeed, although these techniques may result in better disentangling of the manifold space, this often comes at the expense of generating complex and less interpretable features that cannot be related to brain tissue microstructure.

Overall, data representation frameworks such as the one presented here can become fundamental in advancing the application of diffusion models in health and disease. The proposed framework may open new avenues for examining brain microstructure in general and other related lines of research, especially if complimented with other modalities such as measures derived from quantitative magnetisation transfer (Rovaris et al., 2003; Cercignani and Bouyagoub, 2018). Indeed, with the ever-growing acquisition of large cohorts of subjects, feature extraction techniques may become essential tools for processing multi-dimensional brain imaging datasets (e.g., the Human Connectome Project with >1000 young adults scans; Van Essen et al. (2013) or the UK Biobank with its 500,000 participants; Miller et al. (2016)).

Finally, our study may also open new avenues for fibre clustering by leveraging microstructural properties mapped over fibre bundles. Suppl. Fig. 1 shows how different bundles project and cluster in the new reference frame formed by PC1 and PC2. One can observe that PC1 is sensitive to various hindrance level in white matter by disentangling bundles such as the CST (green), genu (blue) and splenium (pink). Conversely, pathways that are known to have many crossing regions such as the AF and the SLF are located on the superior portion of the bi-plot, showing properties of increased orientational dispersion (Suppl. Fig. 1 orange and purple, respectively).

Conclusions

In summary, our findings demonstrate that there exist redundancies in measures conventionally derived from dMRI and that those redundancies may be exploited by dimensionality reduction to reduce the risk of Type I errors, arising from multiple statistical comparisons. Our results support the use of data reduction to detect along-tract differences in tissue microstructure. Specifically, the curse of dimensionality and redundancies in statistical analyses were considerably mitigated by extracting principal components that summarise the inter-dependent measures. From an application perspective, a general increase in the first component related to white matter hindrance was found to have a significant correlation with age in various developmentally sensitive pathways, a change that would otherwise remain undetected using conventional approaches.

Funding

MC is supported by the Postdoctoral Fellowships Program from the Natural Sciences and Engineering Research Council of Canada (PDF-502385-2017) and a Wellcome Trust New Investigator Award (to DKJ). ER is supported by a Marshall Sherfield Postdoctoral Fellowship. CMWT is supported by a Rubicon grant from the Netherlands Organisation for Scientific Research (680-50-1527). This work was also supported by a Wellcome Trust Investigator Award (096646/Z/11/Z), a Wellcome Trust Strategic Award (104943/Z/14/Z), and an EPSRC equipment grant (EP/M029778/1).