Computing and visualising intra‐voxel orientation‐specific relaxation–diffusion features in the human brain

Abstract Diffusion MRI techniques are used widely to study the characteristics of the human brain connectome in vivo. However, to resolve and characterise white matter (WM) fibres in heterogeneous MRI voxels remains a challenging problem typically approached with signal models that rely on prior information and constraints. We have recently introduced a 5D relaxation–diffusion correlation framework wherein multidimensional diffusion encoding strategies are used to acquire data at multiple echo‐times to increase the amount of information encoded into the signal and ease the constraints needed for signal inversion. Nonparametric Monte Carlo inversion of the resulting datasets yields 5D relaxation–diffusion distributions where contributions from different sub‐voxel tissue environments are separated with minimal assumptions on their microscopic properties. Here, we build on the 5D correlation approach to derive fibre‐specific metrics that can be mapped throughout the imaged brain volume. Distribution components ascribed to fibrous tissues are resolved, and subsequently mapped to a dense mesh of overlapping orientation bins to define a smooth orientation distribution function (ODF). Moreover, relaxation and diffusion measures are correlated to each independent ODF coordinate, thereby allowing the estimation of orientation‐specific relaxation rates and diffusivities. The proposed method is tested on a healthy volunteer, where the estimated ODFs were observed to capture major WM tracts, resolve fibre crossings, and, more importantly, inform on the relaxation and diffusion features along with distinct fibre bundles. If combined with fibre‐tracking algorithms, the methodology presented in this work has potential for increasing the depth of characterisation of microstructural properties along individual WM pathways.


S.1 Diffusion-encoding gradient waveforms
Figure S1.Set of gradient waveforms used in this study.The ST acronym identifies a standard Stejskal-Tanner waveform whose spectral profile (Callaghan & Stepišnik, 1996;Lundell et al., 2019) is distinct from those of the non-monopolar waveforms.The waveforms yielding b D = -0.5, 0, and 0.5 were optimized according to the numerical procedure described in refs.(Sjölund et al., 2015) and (Szczepankiewicz, Westin, & Nilsson, 2019).The displayed waveforms were inserted within a spinecho sequence with an EPI readout.To clarify the locations of the spin-echo radio-frequency pulses and the EPI block, we divide each waveform in two components: G 1 (t) and G 2 (t).The 90° pulse is executed at t = 0, before the G 1 (t) component, while the 180° pulse is applied after a time t = t E /2 and is bracketed by the G 1 (t) and G 2 (t) components.Signal readout starts shortly after the conclusion of G 2 (t), with k = 0 of the EPI block coinciding with the echo-time t E .Relaxation-weighting is achieved by varying t E while keeping constant the location of G 1 (t) and G 2 (t) relative to the 180° pulse.

S.2 Orientation dependence of R 2 -D peak-specific metrics
To explore a possible angular dependence of the R 2 -D metrics from different fibre populations, we computed peak-specific Ê[X] values and compared them against their respective q coordinates.The results are displayed in Figure S2, where no clear relationship between Ê[X] and peak orientation is observed for any of the extracted metrics.This observation is in contrast with previous in vivo brain MRI studies where a relationship between the R 2 rates of myelinated fibres and their orientation relative to B 0 (Gil et al., 2016;Knight, Wood, Couthard, & Kauppinen, 2015;McKinnon & Jensen, 2019;Tax, Kleban, Barakovic, Chamberland, & Jones, 2020) has been reported.In particular, Gil and co-workers (Gil et al., 2016) have estimated an angular variation of DR 2 ~ 1.5 s -1 for healthy white matter (WM) tissue.Figure 2A of the main document shows that the uncertainty of the Monte Carlo analysis procedure can introduce R 2 differences of up to DR 2 ~ 2.4 s -1 within a single fibre population; this suggests that a subtle R 2 variation is very challenging to resolve with our minimally constrained approach and explains the approximately constant value observed in Figure S2 for the Ê[R 2 ] and Ê[T 2 ] metrics.To better assess the minimum R 2 uncertainty of our analysis protocol, we focused on individual CC voxels yielding single-lobe ODFs and computed the interquartile range of mean R 2 values obtained from different bootstrap ℇ  b thin components.This yielded interquartile ranges of 1.5-2.5 s -1 for the various CC voxels, providing further evidence that the R 2 uncertainty of a single fibre population determined through the Monte Carlo inversion is on the same order of magnitude as the R 2 orientational effects reported in ref. (Gil et al., 2016).

Figure S2.
Peak-specific means, Ê[X], of T 2 , R2, isotropic diffusivity D iso , and squared normalized diffusion anisotropy D D 2 plotted as a function of q, the polar angle describing the orientation of the various peaks relative to the laboratory frame of reference.The peaks were estimated from local maxima of the smooth ODF, as described in the Methods section.The q angles were sorted into 30 different bins.The solid grey, solid black, and dashed grey lines represent the 75 th , 50 th (or median), and 25 th percentile, respectively, of each angle bin.The shaded grey region illustrates the interquartile range of each angle bin.

S.3 R 2 -D coloured tractography
To explore the potential of the suggested framework in teasing out WM pathways with distinct R 2 -D properties, we implemented a novel visualization approach in FiberNavigator (Chamberland, Whittingstall, Fortin, Mathieu, & Descoteaux, 2014) allowing for the colouring of tracks according to their respective peak amplitudes during real-time tracking.T 2 and D D 2 values were subsequently encoded in the amplitudes of the various peaks, and the new FiberNavigator feature used to map the the derived tracks.The results of this approach are displayed in FigureS3, which shows a colour-coded tractogram where T 2 and D D 2 values are mapped along each streamline.Future, improved versions of the tractograms displayed in FigureS3may be used to investigate differences in the microstructural properties of crossing WM pathways.

Figure S3 .Figure S4 .
Figure S3.Tractogram with streamlines coloured according to orientation-resolved means, Ê[X], of T2 (left panel) and squared normalized diffusion anisotropy D D 2 (right panel).Top: axial view; Middle: coronal view; Bottom: sagittal view.Global opacity rendering was applied to improve visualisation and remove occlusion.