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.

Keywords: diffusion MRI; fibre ODF; fibre-specific metrics; partial volume effects; tensor-valued diffusion encoding; white matter.

FIGURE 1

Resolution of sub‐voxel fibre‐like components and subsequent estimation of the associated colour‐coded Orientation Distributions Functions (ODFs). (a) R ₂‐D distribution obtained for a voxel containing both CSF and two crossing WM populations. The 5D P(R ₂,D) is reported as a 3D logarithmic scatter plot of R ₂, isotropic diffusivities D _iso, and axial–radial diffusivity ratios D _‖/D _⊥, with circle area proportional to the weight of the corresponding R ₂‐D component, w. Colour coding is defined as: [R,G,B] = [cosφ sinθ, sinϕ sinθ, cosθ] · |D _‖−D _⊥/max(D,D _⊥), where (θ,ϕ) gives the orientation of each axisymmetric D. The R ₂‐D space is divided into three coarse bins named ‘big’ (blue volume), ‘thin’ (red volume), and ‘thick’ (green volume). Components falling in the ‘thin’ bin are singled‐out and interpreted as fibres. (b) Spatial distribution of per‐bin signal contributions. The middle map shows the fractional populations in the ‘big’ (blue), ‘thin’ (red), and ‘thick’ (green) bins as a colour‐coded composite image. The rightmost map focuses on the signal contributions from components within the ‘thin’ subset, f _thin, the complement of which, (1 − f _thin), gives the signal fraction from all components not used for ODF calculation. The crosses locate the voxel whose distribution is shown in panel (a). (c) Scheme for calculating colour‐coded ODFs. The R ₂‐coloured circles denote the ‘thin’ components from a bootstrap solution of the voxel signalled in panel (b). Circle area is proportional to w, while the [x,y,z] circle coordinates are defined as either [cosϕ sinθ, sinϕ sinθ, cosθ] (left) or [cosϕ sinθ, sinϕ sinθ, cosϕ] · w (middle and right). In the left plot, the discrete R ₂‐D components are displayed on a unit sphere represented by a 1,000‐point (θ,ϕ) mesh. The weights of the P(R ₂,D) components are first mapped to the mesh through Equation (6), creating an ODF glyph whose radii inform on the R ₂‐D probability density along a given (θ,ϕ) direction (middle). Following the ODF estimation, Equation (9) is used to assign mean values of R ₂,D _iso, or D _Δ to each mesh point and define the colour the ODF glyph (right)

FIGURE 2

R ₂‐D distributions and Orientation Distribution Functions (ODF) retrieved for in silico fibre‐crossing datasets. (a) 5D P(R ₂,D) distributions displayed as 2D scatter plots of log(R ₂) and θ, the polar angle defining D orientation. Circle area is proportional to the weight of the corresponding component and colouring is defined as R,G,B] = [cosϕ sinθ, sinϕ sinθ, cosθ] · |D _‖ − D _⊥/max(D _‖,D _⊥), where D _‖ and D _⊥ denote the axial and radial diffusivities, respectively, and ϕ is the azimuth angle of D. The yellow crosses identify the ground‐truth values. (b) ODF glyphs estimated from the distributions in panel (a), using Watson kernel with different orientation dispersion factors (see Equation (6) for further details). The ODF colouring follows a conventional directional scheme: [R,G,B] = [μ _xx,μ _yy,μ _zz], where μ _ii are the elements of the unit vector μ(θ,ϕ) defining the orientation of mesh‐point (θ,ϕ). (c) ODF glyphs estimated for an in vivo voxel rendered as triangular surface plots for varying angular standard deviation σ of the convolution kernel at constant number of triangle vertices and the same underlying 5D P(R ₂,D) distribution. The depicted voxel comprises a two‐way crossing between fibres from the corticospinal tract and the corpus callosum

FIGURE 3

Orientation Distribution Functions (ODF) and peak metrics retrieved for in silico fibre‐crossing datasets. (a) ODF glyphs computed for the in silico systems described in Section 2.6. The displayed ODFs are coloured according to the orientation‐resolved means of T ₂, Ê[T ₂] (see Equation (9) for more details). The inset displays the f _thin = 0.3 ODFs with their base amplitude rescaled by a factor of three. (b–d) Number of peaks (# peaks) and mean biases (E[ΔX]) in peak angle (b), peak T ₂, and peak squared normalised diffusion anisotropy $D_{\mathrm{\Delta }}^2$ estimated for three selected in silico systems. We refer the reader to Section 2.6 for a detailed description of the designed systems. Briefly, the displayed systems were designed by combining upto three discrete fibre components (D _iso = 0.75 × 10^–9 m² s^–1, D _Δ = 0.9) with different signal fractions (1 − f _thin) of a GM‐like component (T ₂ = 90 ms, D _iso = 0.8 × 10^–9 m² s^–1, D _Δ = 0.2). Here, we inspect a set of n _fibre discrete fibres at various crossing angles θ _cross (b), two‐fibre crossings for various T ₂ values of the crossing fibre, T _2,cross (c), and two‐fibre crossings with varying f _thin (d). In (b) and (c), we set 1 − f _thin = 0.3, while in (b) and (d), we selected T _2,cross = 80 ms. The metrics were calculated across 40 different noise realisations; the points represent the mean over the various signal realisations, while the error bars indicate the standard deviation across signal realisations. The dashed grey line indicates the true number of peaks (left plots) or the zero‐bias line

FIGURE 4

Per‐voxel Orientation Distribution Functions (ODF), P(θ,ϕ), estimated from R ₂‐D distribution components ascribed to the ‘thin’ bin defined in Figure 1. The voxel‐wise P(θ,ϕ) were computed by using Equation (6) to map the weights of the bin‐resolved discrete P(R ₂,D) components into a 1,000‐point spherical mesh. Here, each ODF is represented as a 3D polar plot with a local radius given by P(θ,ϕ) and colour‐coded according to [R,G,B] = [μ _xx,μ _yy,μ _zz], where μ _ii are the elements of the unit vector μ(θ,ϕ) (see Equation (6) for further details). In the left and top‐right panels, the sets of ODF glyphs are superimposed on a grey‐scaled map that shows the signal contributions from non‐fibre‐like components (1 − f _thin), that is, signal fractions from the ‘big’ and ‘thick’ populations. The zoom‐ins in the lower‐right panel offer a more detailed look into selected fibre crossing regions (continuous line boxes) and three‐fibre crossing voxels (dashed line boxes) found in the centrum semiovale. The various arrows identify fibre tracts mentioned in the main text

FIGURE 5

Opacity rendering of streamline tractography data, where the opacity reflects streamline density (computed using a slice thickness of 3 voxels). The various tracks are coloured according to their orientation: red (left‐right), green (anteroposterior), and blue (superoinferior). The right‐side panels display one coronal slice (top) and two different sagittal slices (middle and bottom), while the left‐side panel displays an axial slice

FIGURE 6

Orientation Distribution Function (ODF) maps coloured according to orientation‐resolved means, Ê[X], of R ₂, isotropic diffusivity D _iso, and squared normalised diffusion anisotropy $D_{\mathrm{\Delta }}^2$. All Ê[X] were calculated using Equation (9) and are displayed on a linear scale. The lower panel displays a zoom into a region containing fibre crossings between the corpus callosum and the corticospinal tract. The dashed‐line box in the top‐left map identifies the high‐R ₂ fibres found in the globus pallidus

FIGURE 7

Orientation Distribution Function (ODF) maps coloured according to the orientation‐resolved means of T ₂, Ê[T ₂]. The Ê[T ₂] values are displayed on a linear colour scale. The left and top‐right panels display the sets of ODF glyphs superimposed on a grey‐scaled map showing the signal fractions from the ‘big’ and ‘thick’ bin populations (1 − f _thin) (non‐fibre‐like components). The zoom‐ins in the lower‐right panel offer a more detailed look into selected regions (continuous line boxes) and voxels (dashed line boxes) containing crossing between fibre populations with distinct Ê[T ₂]. The observed high‐T ₂ components are assigned to the forceps major (yellow boxes) and the corticospinal tract (magenta boxes)

FIGURE 8

Orientation‐resolved metrics estimated for a two‐fibre‐crossing voxel in the superior longitudinal fasciculus. (a) Orientation Distribution Function (ODF) estimated for the selected voxel. The black points identify the two peaks of the displayed ODF, peaks A and B. (b,c) Fibre‐specific R ₂‐D metrics. The (θ,ϕ) orientation space was divided into four quadrants centred on A, B, and their corresponding antipodes; ‘thin’ R ₂‐D components, ${\xi }_{n_{\mathrm{b}}}^{\text{thin}}$, were then assigned to either fibre population A or fibre population B depending on their (θ,ϕ) coordinates (e.g., ${\xi }_{n_{\mathrm{b}}}^{\text{thin}}$ components falling into the quadrant centred on peak A, are assigned to fibre population A). For each orientation bin and each bootstrap, we estimate the mean signal fraction, R ₂, isotropic diffusivity D _iso, squared normalised diffusion anisotropy $D_{\mathrm{\Delta }}^2$, and orientation, thus obtaining a set of 96 × 6 scalars: 96 different estimates of six distinct parameters. (b) Ensemble of fibre‐resolved orientations displayed on the unit sphere. The colouring of the sphere identifies the (θ,ϕ) space assigned to each fibre population. The coloured lines indicate the peak orientation of fibres A (green) and B (red), while the black lines indicate the [x,y,z] coordinates. (c) Boxplots displaying the average and dispersion of the fibre‐resolved signal fractions, R ₂,D _iso, and $D_{\mathrm{\Delta }}^2$. The average was estimated as the median, while dispersion was assessed as the interquartile range. The whiskers identify the maximum and minimum estimated values