Visualizing orientation-specific relaxation-diffusion features mapped onto orientation distribution functions estimated via nonparametric Monte Carlo MRI signal inversion

Diffusion MRI techniques are widely used to study in vivo changes in the human brain connectome. However, to resolve and characterise white matter fibres in heterogeneous diffusion MRI voxels remains a challenging problem typically approached with signal models that rely on prior information and restrictive 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 in order 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 in order 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 distinct fibre bundles. If combined with fibre-tracking algorithms, the methodology presented in this work may be useful for investigating the microstructural properties along individual white matter pathways.


MAIN TEXT INTRODUCTION
For decades, neuroscience has focused not only on unravelling the function of brain areas (Brodmann, 1909;Zilles & Amunts, 2010), but also on the communication between these areas (Sporns, Tononi, & Kötter, 2005). Brain function fundamentally depends on the effective transport of information (Passingham, Stephan, & Kötter, 2002), and inter-individual functional differences may well arise from differences in composition and structure of the physical connections (Alexander-Bloch, Giedd, & Bullmore, 2013). Studying the constituents of the white matter (WM) that form these connections is therefore of paramount importance to understanding brain function in health, disease and development. The advent of diffusion MRI techniques, which can probe structures at much smaller scales than the imaging resolution by virtue of sensing the random motion of water molecules, has undoubtedly increased the interest in studying WM in the living brain. The possibility to derive quantitative features sensitive to tissue microstructure Le Bihan, 1995), and to virtually reconstruct brain connections with fibre-tracking algorithms (Basser, Pajevic, Pierpaoli, Duda, & Aldroubi, 2000;Mori, Crain, Chacko, & Van Zijl, 1999) led to the quick adoption of diffusion MRI in many clinical research applications (Barnea-Goraly et al., 2004;Lebel, Walker, Leemans, Phillips, & Beaulieu, 2008;Lim et al., 1999;Werring, Clark, Barker, Thompson, & Miller, 1999). More recently, tractometry techniques have been developed to tease out WM pathways and characterize their individual tissue microstructure by mapping sets of diffusion-derived parameters along the extracted tracks (Bells et al., 2011;Chamberland et al., 2019;De Santis, Drakesmith, Bells, Assaf, & Jones, 2014;Rheault, Houde, & Descoteaux, 2017;Yeatman, Dougherty, Myall, Wandell, & Feldman, 2012). Fibre-tracking techniques typically rely on the estimation of a fibre Orientation Distribution Function (ODF) per voxel, which is a function on the unit sphere representing the fraction of fibres in each direction (Dell'Acqua & Tournier, 2019;Tournier, 2019). It should be noted that the influence of the fibre ODF is distinct from the orientation distribution of the diffusion signal, and its extraction relies on assessing how tissue microstructure influences the measured MRI signal. Diffusion MRI studies of WM commonly assume that the voxel-level microstructural features can be adequately represented by a single signal response function (Dell'Acqua & Tournier, 2019;Novikov, Fieremans, Jespersen, & Kiselev, 2019). Under this assumption, the measured signal is written as the convolution between the fibre ODF and a kernel describing the signal response of a set of fibres with a common orientation. The simultaneous unconstrained estimation of the ODF and the microstructural kernel, however, has shown to be notoriously challenging for the diffusion MRI protocols typically used for in vivo research studies (Jelescu, Veraart, Fieremans, & Novikov, 2016). The complexity of this problem is commonly reduced by imposing a set of priors and constraints. Spherical deconvolution of the diffusion MRI signal (Anderson, 2005;Dell'Acqua et al., 2007;Dell'Acqua & Tournier, 2019;Tournier, Calamante, & Connelly, 2007;Tournier, Calamante, Gadian, & Connelly, 2004), for example, determines an empirical kernel for the whole brain representing the signal response of a single-fibre population and subsequently solves for the ODF. For heterogeneous voxels containing not only WM but also unknown amounts of grey matter (GM), cerebrospinal fluid (CSF), or pathological tissue, this approach can yield biased ODF estimates. Multi-tissue spherical deconvolution (Jeurissen, Tournier, Dhollander, Connelly, & Sijbers, 2014) has been proposed to simultaneously resolve sub-voxel tissue fractions and the fibre ODF. While this technique can be used to separate the sub-voxel signal contributions from WM, GM, and CSF, it still assumes a single kernel for all voxels of a given tissue type, which needs to be calibrated a priori (Tax, Jeurissen, Vos, Viergever, & Leemans, 2014). Inaccuracies of the calibrated kernels can further bias the estimated fractions and fibre ODFs (Guo et al., 2019;Parker et al., 2013). Alternatively, the voxel-wise kernel can be estimated by first factoring out the ODF through the computation of rotational invariants, and then fitting the data to signal models that set a pre-defined number of microscopic environments with potentially constrained diffusion properties (Kaden, Kelm, Carson, Does, & Alexander, 2016;Novikov et al., 2019;Novikov, Veraart, Jelescu, & Fieremans, 2018). However, different fibre populations within a voxel likely exhibit different microstructural properties (Aboitiz, Scheibel, Fisher, & Zaidel, 1992;De Santis, Assaf, Jeurissen, Jones, & Roebroeck, 2016;Howard et al., 2019;Scherrer et al., 2016), which cannot be reflected with a single voxel-wise kernel. It should furthermore be noted that transverse relaxation time T2 differences between distinct tissue types are often ignored, which can further bias the quantification of tissue fractions with a single fibre response kernel. The possible existence of a variation of T 2 in anisotropic structures with respect to the orientation of the main magnetic field B 0 (Lindblom, Wennerström, & Arvidson, 1977), well known in studies of cartilage structure (Henkelman, Stanisz, Kim, & Bronskill, 1994) and more recently reported in in vivo human WM studies (Gil et al., 2016;Knight, Wood, Couthard, & Kauppinen, 2015;McKinnon & Jensen, 2019), would introduce an additional T 2 dispersion and further complicate the quantification of sub-voxel signal fractions. The possible existence of T 2 differences between distinct fibre bundles has motivated the recent development of methods allowing for the measurement of fibre-specific estimates of the transverse relaxation time (de Almeida Martins & Topgaard, 2018;Ning, Gagoski, Szczepankiewicz, Westin, & Rathi, 2020;Schiavi et al., 2019). Inspired by multidimensional solid-state NMR methodology (Schmidt-Rohr & Spiess, 1994;Topgaard, 2017), we have introduced a framework to quantify the composition of each voxel with joint distributions of effective transverse relaxation rates R2 = 1/T2 and apparent diffusion tensors D de Almeida Martins & Topgaard, 2018). Specifically, the inclusion of diffusion MRI data measured with multidimensional diffusion encoding schemes (Topgaard, 2017) and different echo times was observed to alleviate the constraints needed to resolve sub-voxel tissue heterogeneity . Capitalizing on these acquisitions, we quantified sub-voxel compositions using 5D discrete R2-D distributions retrieved from the data using a nonparametric Monte Carlo inversion procedure. However, visualising the retrieved sub-voxel information is challenging because of the high dimensionality of the distributions. The challenge of visualizing the intricate and comprehensive information within diffusion MRI datasets is an active area of research (Leemans, 2010;Schultz & Vilanova, 2019) and very well established visualization strategies exist to either convey the tensorial properties of a single voxel-averaged D (Kindlmann, 2004;Westin et al., 1999) or to visualize a continuous ODF (Peeters, Prckovska, Almsick, Vilanova, & Romeny, 2009;Schultz & Kindlmann, 2010;Tournier et al., 2004;Tuch et al., 2002). However, such techniques are not immediately applicable to the discrete multi-component distributions retrieved with our 5D correlation framework. Previously, we converted the retrieved distributions to sets of statistical parameter maps derived from either the entirety or sub-divisions ('bins') of the distribution space Topgaard, 2019). In ref. , the R2-D space was divided into three bins capturing different ranges of D eigenvalues in order to separate the signal contributions from microscopic tissue environments with distinct diffusion properties. Even though bin-resolved maps of signal fractions and means were observed to be useful to map sub-voxel heterogeneity throughout the imaged brain volume, they do not provide information on orientation-resolved properties. In this contribution, we demonstrate how R2-D distributions can be used to derive and visualize fibre-specific relaxation and diffusion metrics. This is done by extending the binning procedure to the space of D orientations, and mapping discrete P(R2,D) components to a spherical mesh representing a dense set of orientation bins. The attained orientation-resolved information is then conveyed as maps of ODF glyphs that are colour-coded according to the underlying relaxation and diffusion properties; this greatly facilitates the inspection and interpretation of the orientational variation of the 5D P(R2,D). The representation of the high-dimensional information as colour-coded ODFs is furthermore compatible with tractography algorithms which hence allows the extension to visualisation of longer-range properties in 3D.

Estimation of 5D relaxation-diffusion distributions
In diffusion MRI, heterogeneous tissues can be described as a collection of microscopic tissue environments wherein water diffusion is modelled by a local apparent diffusion tensor D. Within this multi-tensor approach, the diffusion MRI signal is approximated as a weighted sum of the signals from the individual microscopic tissue environments (Jian, Vemuri, Özarslan, Carney, & Mareci, 2007;Novikov et al., 2019;Westin et al., 2016). A similar heterogeneous description has also been used in R 2 studies of intra-voxel brain tissue structure (Does, 2018;MacKay et al., 2006). The transverse relaxation signal of water within tissues is typically expressed as a multi-exponential decay, given by the Laplace transform of a probability distribution of R 2 values (Kroeker & Mark Henkelman, 1986;Whittall & MacKay, 1989;Whittall et al., 1997). Each coordinate of the relaxation probability distribution characterizes the signal fraction of the microscopic environment with the corresponding R 2 rate. Combining the relaxation and diffusion descriptions, the detected signal S(t E ,b) can be written as where P(R 2 ,D) is the continuous joint probability distribution of R 2 and D, t E denotes the echo-time, b is the diffusion-encoding tensor, and S 0 is the signal amplitude at vanishing relaxation-and diffusion-weighting, i.e. S 0 = S(t E = 0, b = 0). The integration of D spans over the space Sym + (3) of symmetric positive semi-definite 3×3 tensors. The kernel K(t E ,b,R 2 ,D) encapsulates the signal decays mapping the distribution onto the detected signal. Constraining the integral in Eq.
(1) to the space of axisymmetric diffusion tensors, each D can be parameterized by its axial and radial diffusivities, D || and D^, and by the polar and azimuthal angles, q and f, that define its orientation. The D || and D^ eigenvalues can in turn be combined to define measures of isotropic diffusivity D iso = (D || + 2D^)/3 and normalized diffusion anisotropy D D = (D || -D^)/(3D iso ) (Eriksson, Lasič, Nilsson, Westin, & Topgaard, 2015). Using the popular approach of approximating the signal from each microscopic environment as an exponential decay (Dell'Acqua et al., 2007;Does, 2018;Kaden et al., 2016;MacKay et al., 2006;Novikov et al., 2018;Scherrer et al., 2016;Tuch et al., 2002;Westin et al., 2016), considering only axisymmetric b, and adopting the (D iso ,D D ,q,f) parametrization, Eq. (1) can be expanded as (de Almeida Martins & Topgaard, 2018) where b = Tr(b) is recognized as the traditional b-value and b D denotes the normalized anisotropy of the diffusion-encoding tensor (Eriksson et al., 2015). P 2 (x) = (3x 2 -1)/2 is the second Legendre polynomial, and b is the shortest angle between the symmetry axes of D and b. Note that each diffusion orientation (q,f) is associated to its own set of microscopic properties (R 2 ,D iso ,D D ) and that no overarching microstructural kernel or universal orientation structure is assumed. This means that Eq.
(2) allows for fibre populations with distinct R 2 -D properties. For numerical implementation, Eq.
(2) is discretized as s = Kw, where s is the column vector of signal amplitudes measured with M combinations of (t E ,b) values, K is the inversion kernel matrix, and w is a vector containing the weights w n of N discrete (R 2,n ,D ||,n ,D^, n ,q n ,f n ) configurations. The estimation of w can then be cast as a constrained linear least-squares problem w = arg min w≥0 ||s − Kw|| 2 2 . (4) In practice, the argument-minimum operator in Eq. (4) is replaced by a softer condition that searches for a solution within the noise variance. While seemingly straightforward, the problem of finding a solution whose residuals are compatible with the experimental noise is a poorly conditioned one. Indeed, multiple distinct solutions can be found to fit a single noisy dataset. This has motivated the development of several regularization strategies in order to improve the stability of the inverse problem (Daducci et al., 2015;Mitchell, Chandrasekera, & Gladden, 2012;Provencher, 1982;Whittall & MacKay, 1989). A common strategy is to incorporate a regularization term that promotes either a smooth (Benjamini & Basser, 2017;Provencher, 1982;Slator et al., 2019;Venkataramanan, Song, & Hurlimann, 2002) or a sparse (Benjamini & Basser, 2016;Berman, Levi, Parmet, Saunders, & Wiesman, 2013;Urbańczyk, Bernin, Koźmiński, & Kazimierczuk, 2013) w solution at the expense of a higher residual error. Monte Carlo algorithms have been used in the porous media field as an alternative to conventional regularized approaches (de Almeida Martins & Topgaard, 2016de Kort, van Duynhoven, Hoeben, Janssen, & Van As, 2014;Prange & Song, 2009). These algorithms purposely explore the variability between solutions and estimate an ensemble of distributions consistent with the experimental data. In this work, we use an unconstrained Monte-Carlo algorithm specially designed to handle high-dimensional correlation datasets (de Almeida Martins & Topgaard, 2018;Topgaard, 2019). The algorithm can be broadly divided in two iteration cycles. In the first cycle, the proliferation cycle, a user-defined N in number of points is randomly selected from the (log(R 2 ),log(D || ),log(D^),cosq,f) space, and the corresponding set of weights is found by solving Eq. (4) via a non-negative linear leastsquares algorithm (Lawson & Hanson, 1974); points with non-zero weights are stored and merged with a newly-generated random set. This procedure is repeated for a total of N p times, and N p random sets of (R 2,n ,D ||,n ,D^, n ,q n ,f n ) components are sampled in order to find a configuration yielding sufficiently low residuals. The resulting {(R 2,n ,D ||,n ,D^, n ,q n ,f n )} configuration is stored, duplicated, and its duplicate is then subjected to a small random perturbation. This initiates the second iteration cycle, named mutation cycle, wherein configurations compete with their perturbed counterparts on the basis of lowest sum of squared residuals. The mutation cycle comprises a number of N m rounds, following which a possible solution is estimated by selecting the points with the N highest weights. In this work we sampled N in = 200 points from the (0 < log(R 2 / s -1 ) < 1.5, -11.3 < log(D || / m 2 s -1 ) < -8.3, -11.3 < log(D^ / m 2 s -1 ) < -8.3, 0 < cosq < 1, 0 < f < 2p) space, and used N p = 20, N m =20 , and N = 20. This inversion was performed voxel-wise and bootstrap with replacement was used in order to estimate per-voxel ensembles of N b = 96 solutions, each of which consisting of 20 (R 2,n ,D ||,n ,D^, n ,q n ,f n ) components, {(R 2,n ,D ||,n ,D^, n ,q n ,f n )} 1£n£N=20 , and their respective w n weights.

Resolution of sub-voxel fibre components
Spatially-resolved 5D R 2 -D distributions were estimated using the procedure described in the previous section. As the main brain components -white matter (WM), grey matter (GM), and cerebrospinal fluid (CSF) -are characterized by clearly distinct diffusion properties, we expect most (R 2,n ,D ||,n ,D^, n ,q n ,f n ) components to agglomerate within three distant regions of the diffusion space (Pierpaoli, Jezzard, Basser, Barnett, & Di Chiro, 1996). The idea that most P(R 2 ,D) components will fall within three coarse regions has inspired the division of the R 2 -D space into three smaller subsets ('bins') based on the diffusion properties of WM, GM, and CSF . We then defined three bins named 'thin' (0.6 < log(D || /D^) < 3.5, -10 < log(D iso /m 2 s -1 ) < -8.7, -0.5 < log(R 2 /s -1 ) < 2), 'thick' (-3.5 < log(D || /D^) < 0.6, -10 < log(D iso /m 2 s -1 ) < -8.7, -0.5 < log(R 2 /s -1 ) < 2), and 'big' (-3.5 < log(D || /D^) < 3.5, -8.7 < log(D iso /m 2 s -1 ) < -8, -0.5 < log(R 2 /s -1 ) < 2). The names of the different bins highlight the geometry of the D captured by each one of them. For each bootstrap realization n b (1 £ n b £ N b ), signal contributions from anisotropic tissues are resolved by selecting the set of P(R 2 ,D) components that fall within the 'thin' bin: as representing the R 2 -D properties and signal fractions, respectively, of a discrete set of sub-voxel fibre populations. The binning and anisotropic selection processes are illustrated in panels A and B of Figure 1.

ODF estimation
The colour-coded 3D scatter plots of R 2 , D iso , and D || /D^ displayed in Figure 1 allow the visualization of the full set of properties of the voxel-wise ℇ b thin components. Despite its usefulness, the scatter plot representation concentrates points corresponding to anisotropic components within a small region of the (R 2 , D iso , D || /D^) space, which in turn makes it difficult to evaluate their orientation properties in detail. For example, while Figure 1A clearly informs on the existence of two fibre populations oriented along two different directions (red and green points), it does not provide unambiguous information about the relative signal contributions of the two populations. To better inspect the orientational information of the underlying P(R 2 ,D), it is helpful to convert the discrete set of fibre orientations to a continuous object informing on the R 2 -D probability density in each direction, which can then be visualized as a single glyph with an intuitive geometrical interpretation. To achieve that, we used a 1000-point triangle mesh on the unit sphere, created via an electrostatic repulsion scheme (Bak & Nielsen, 1997;Jones, Horsfield, & Simmons, 1999), to create 1000 orientation bins; the mesh vertices define the centres of (q,f) bins without clearly defined boundaries. The median angular distance between nearest-neighbouring mesh points (or bin centres) is approximately 7°. Afterwards, a smoothing kernel was used to map the weights of ℇ b thin onto the dense set of (q,f) bins. The role of the smoothing kernel is to weight the influence of each bin according to the angular distance between its centre and a given ℇ b thin component, and to distribute the contributions of each discrete {w i } iÎ{n b ,thin} throughout various bins in order to define a smooth Orientation Distribution Function (ODF) P n b (q,f) that can be straightforwardly visualized as a polar plot (Leemans, 2010;Schultz & Vilanova, 2019). In this work, the orientation binning and consequent estimation of P n b (q,f) was performed through a convolution with a Watson distribution kernel (Mardia & Jupp, 2009;Watson, 1965): where u i is the unit vector describing the orientation of the i-th discrete component and µ (q,f) is the unit mean direction vector of bin centre (q,f). The variable k denotes a concentration parameter that regulates the amount of dispersion about µ (q,f). Further insight onto the nature of the Watson distribution kernel and the role of parameter k is attained by considering the small-angle approximation of the former where b is the shortest angle between unit vectors u and µ. Within this approximation, the Watson kernel is rewritten as a familiar Gaussian smoothing kernel whose standard deviation, s, defines an angular spreading that is directly related to the concentration parameter k, s = (2k) -1/2 . The relationship between k and s enables us to easily set a dispersion factor in relation to the angular distance of the various mesh points. To cover the minimum angular distance between mesh points, we set s to be 50% higher than the median angular distance between nearest-neighbouring mesh points, i.e. s = 10.5° or, equivalently, k = 14.9. The rationale behind our choice of k is elaborated upon in section 3.1. Separate ODFs were calculated for each of the N b = 96 bootstrapped P(R 2 ,D) solutions; the final P(q,f) was then estimated as the median of the N b independent ODFs: As the orientation of each fibre solution is correlated to a given set of {R 2,n ,D ||,n ,D^, n } coordinates, we can assign any statistical descriptor of the (R 2 ,D iso ,D Δ ) space to the various coordinates of P(q,f).
In line with previous works Topgaard, 2019) where R 2 -D were converted into bin-resolved mean maps of R 2 , D iso , and D Δ 2 , we map mean values of R 2 , D iso , and D Δ 2 into the ODF mesh. The mean value of X = T 2 , R 2 , D iso , or D D 2 per mesh orientation for each bootstrap By mapping different descriptors of the R 2 -D distributions to specific ODF coordinates, we can visualize orientation-specific information on tissue composition and structure. As before, the final voxelwise Ê[X](q,f) is estimated as the median of the individual per-bootstrap Ê n b [X](q,f) values: For compactness, we omit the explicit (q,f) dependence from the orientation-resolved means and simply denote them as Ê[X]. The Ê[D Δ 2 ] metric provides orientation-resolved information on the underlying mean-diffusivity. The Ê[D Δ 2 ] metric is the orientation-resolved counterpart of the mean D Δ 2 descriptor (Topgaard, 2019), which is in turn similar to previously introduced anisotropy measures such as the microscopic anisotropy index (MA) (Lawrenz, Koch, & Finsterbusch, 2010), the fractional eccentricity (FE) (Jespersen, Lundell, Sønderby, & Dyrby, 2013), and the microscopic fractional anisotropy (μFA) (Lasič, Szczepankiewicz, Eriksson, Nilsson, & Topgaard, 2014). The mapping of P(R 2 ,D) components to a dense mesh as described by equations (6)-(10) is a key result from this contribution, and provides the basis for extracting and visualizing orientation-resolved information from nonparametric R 2 -D distributions. Figure 1C illustrates how both Ê[X] and the associated orientation distribution can be conveniently represented by colour-coded ODF glyphs; the shape of the glyph reflects the P(q,f) distribution, while the colour informs on the Ê[X] values at the various (q,f) points. Functions used to compute the colour-coded ODFs have been incorporated in the multidimensional diffusion MRI toolbox (Nilsson et al., 2018): https://github.com/JoaoPdAMartins/md-dmri. In this work, maps of ODF glyphs were computed using those same functions and rendered with POV-Ray (http://www.povray.org/). The local maxima of the computed ODFs can be used to estimate "peaks" providing information on the main orientations of the sub-voxel diffusion pattern. The peaks provide a rough and easily accessible measure of the orientation and R 2 -D properties of different ODF lobes and can consequently be used to assess the properties of the various fibre populations. Here, up to four peaks per voxel were determined by assessing the mesh points (q,f) for which P(q,f) is a local maximum and P(q,f)/max (q,f) [P(q,f)] ³ 0.1. The R 2 -D properties of the ODF peaks were estimated by calculating Ê[X] (see Eq. (9)) for each peak orientation.

2.4
In vivo data acquisition Multidimensional relaxation-diffusion MRI data were acquired using a prototype spin-echo diffusion weighted sequence with an echo-planar imaging (EPI) readout, customized for diffusion encoding with user-defined gradient waveforms (Lasič et al., 2014;Szczepankiewicz, Sjolund, Stahlberg, Latt, & Nilsson, 2019). Images were recorded with the following parameters: TR = 4 s, FOV = 234×234×60 mm 3 , voxel-size=3×3´3 mm 3 , partial-Fourier = 6/8, and a parallel-imaging (GRAPPA) factor of 2. Diffusion encoding was performed with a set of five gradient waveforms yielding btensors with four distinct "shapes" (b D = -0.5, 0.0, 0.5, and 1) (Eriksson et al., 2015). The various waveforms were used to probe b-tensors of varying magnitude b, anisotropy b D , and orientation (q,f) at different echo-times t E ; this procedure yields 5D relaxation-diffusion correlated datasets whose dimensions match those of the sought-for distributions. Readers interested in the sequence used in this work are directed to a public repository: https://github.com/filip-szczepankiewicz/fwf_seq_resources. Figure 2 displays the waveforms used in this work and Table 1 summarizes the (t E ,b) acquisition protocol. Besides the (t E ,b) points detailed in Table 1, we additionally acquired b = 0 images with reversed phase-encoding blips at t E = 60, 80, 110, and 150 ms in order to correct for susceptibilityinduced distortions (Andersson, Skare, & Ashburner, 2003). The acquisition scheme comprised a total of 686 (t E ,b) acquired over an imaging session of ~45 mins. The asymmetric waveforms giving b D = -0.5, 0.0, and 0.5 were calculated with a MATLAB package (https://github.com/jsjol/NOW) that optimizes for maximum b (Sjölund et al., 2015) and minimizes the effects of concomitant magnetic field gradients . Linearly encoded (b D = 1) data was acquired with two different gradient waveforms: a non-monopolar gradient waveform and a standard Stejskal-Tanner waveform (Stejskal & Tanner, 1965). The asymmetric gradient pulses from the non-monopolar waveform were manually designed with the aim of minimizing the differences between the spectral profile of that b D = 1 waveform and the spectral profile of the b D ¹ 1 waveforms . The Stejskal-Tanner design was used to probe a shorter t E and higher b-values (b = 4·10 9 m -2 s) than those achievable with the non-monopolar b D = 1 waveform. While the measured apparent diffusivities are known to be related to the frequency spectra of the gradient waveforms (Callaghan & Stepišnik, 1996;Lundell et al., 2019;Stepišnik, 1993), such a relationship is likely to have a negligible effect on healthy human brain data acquired with the limited range of frequency contents probed in this work (Szczepankiewicz, Lasič, et al., 2019) and no biases are expected to originate from the spectral differences of the b D = 1 waveforms. The protocol described above was implemented on a 3T Siemens MAGNETOM Prisma scanner (Siemens Healthcare, Erlangen, Germany) and used to scan a healthy adult volunteer. This study was approved by the Cardiff University School of Psychology ethics committee, and informed written consent was obtained prior to scanning.   (Bak & Nielsen, 1997;Jones et al., 1999) b Repeated for all t E values c Repeated for six different permutations of the [G x ,G y ,G z ] components of the b D = 0 waveform

Post processing
The entire dataset was divided in t E -specific data subsets, which were denoised using random matrix theory , and corrected for Gibbs ringing artefacts using the method described in (Kellner, Dhital, Kiselev, & Reisert, 2016). Signal drift correction was subsequently performed as detailed in ref. (Vos et al., 2017). The acquired data were further corrected for subject motion and eddy-current artefacts using ElastiX (Klein, Staring, Murphy, Viergever, & Pluim, 2009) with extrapolated references (Nilsson, Szczepankiewicz, van Westen, & Hansson, 2015) as implemented in the multidimensional diffusion MRI toolbox (Nilsson et al., 2018); this procedure was performed with the default settings of the toolbox to the entire (t E ,b) dataset. Susceptibility-induced geometrical distortions were corrected using the TOPUP tool in the FMRIB software library (FSL) (Smith et al., 2004), with the same settings being applied to the entire (t E ,b) dataset.

2.6
In silico datasets The angular resolution of our acquisition and analysis protocols was investigated using in silico data. We simulated a two-component system designed to mimic two fibres with similar diffusion features but distinct R2 rates: • Component 1: T 2 = 1/R 2 = 60 ms, D iso = 0.75·10 -9 m 2 s -1 , D D = 0.9, w = 0.5; The ground-truth signal data for the four fibre-crossing systems were generated using the (t E ,b) acquisition scheme indicated in Table 1 and computed from Equation (2). Gaussian distributed noise with an amplitude of 1/SNR was added to the ground-truth signals in order to simulate the effects of experimental noise. The experimental SNR, computed as the mean-to-standard-deviation-ratio of b D = 0 data acquired at b = 0.3·10 9 m -2 s and t E = 80 ms (Szczepankiewicz, Sjolund, et al., 2019), was estimated to SNR = 72 ± 28 for WM regions. Consequently, we defined SNR = 70 for the in silico calculations, a value that is compatible with the SNR of the in vivo data. In line with a recent in silico study of the performance of the Monte Carlo algorithm in inverting multidimensional diffusion data (Reymbaut, Mezzani, de Almeida Martins, & Topgaard, 2020), we drew 100 independent noise configurations and computed 100 different signal realizations for each of the four fibre-crossing systems. The various signal realizations were inverted using the Monte-Carlo algorithm described in section 2.1, and the resulting solution ensembles were subsequently compared against the corresponding ground-truth systems.

Defining the dispersion factor of the Watson kernel
As mentioned in section 2.3, the use of a Watson kernel introduces an artificial angular dispersion to the inverted ℇ b thin components. The amount of angular dispersion is regulated by the user-defined parameter k, and is introduced to assure that the discrete {w i } iÎ{n b ,thin} weights do not fade between the intervals of the µ (q,f) mesh. To understand the smoothing effects of the Watson kernel over a discrete mesh, it is instructive to consider the decay of the Watson function over a given angular distance Db : . Considering a 1000-point mesh and k = 14.9, the values used for the in vivo data analysis, the maximum distance between an arbitrary {q i ,f i } iÎ{n b ,thin} configuration and the nearest mesh point is ~3.5°, a value for which the Watson kernel retains n = 0.95 of its maximal influence. The minimal decay of the Watson kernel over Db = 3.5° ensures that the set of ℇ b thin discrete components is indeed mapped into the mesh. From Equation (7) it is additionally obvious that the choice of k is a trade-off between a sufficiently smooth ODF representation and the angular resolution of the ODF in disentangling different peaks. The question then arises whether setting k = 14.9 may over-smooth the orientational information within the R 2 -D distributions. For instance, with k = 14.9 the Watson kernel will retain more than 50% (n = 0.64) of its maximum value over a distance of Db = 10°, meaning that 20° crossings cannot be resolved with our settings. To assess if the amount of k -generated dispersion is sufficiently low not to misrepresent the orientational information of the R 2 -D distributions, we investigated in silico the angular resolution of the Monte Carlo analysis. The angular resolution of our framework was assessed by inverting in silico data from two anisotropic components crossing at various angles (see section 2.6 for further details). The (R 2 ,q) projections of the attained R 2 -D distributions are displayed in Figure 3A and inform that, at the SNR of the in vivo data, crossings of 30° or higher can be directly resolved in the Monte Carlo P(R 2 ,D). Setting equal R 2 (T 2 = 1/R 2 = 60 ms) properties for both anisotropic components or adding a third isotropic component (T 2 = 1/R 2 = 500 ms, D iso = 2·10 -9 m 2 s -1 , D D = 0) with a total signal fraction of up to 0.2 did not affect the angular resolution of the 30° crossing, but lead to an underestimation of the signal fraction from the q = 30° fibre population. The accurate resolution of the 35° and 40° systems was unaffected by changes in component R 2 or the introduction of the fast-diffusing component. The in silico results then suggest that the maximum achievable angular resolution of our experimental protocol is between 30° and 35°, and a conservative approach is to set k so that 35° crossings are not over-smoothed and obscured. Computing ODFs for the in silico distributions confirms that setting k = 14.9 is indeed sufficient to resolve a 35° crossing (see Figure 3B). While there is room to increase k without risking ℇ b thin disappearing though the holes of the mesh, we observe that a significantly sharper Watson kernel leads to narrow ODF lobes that do not accurately portray the angular dispersion of the underlying R 2 -D distributions (compare panels A and B of Figure 3). Moreover, significantly higher k values were tested in the in vivo dataset and observed to lead to non-smooth 'spiky' ODFs in voxels containing orientationally dispersed fibres.

In vivo fibre orientations
Previous work from our group  has shown that pure component voxels containing either WM, GM, or CSF give rise to clearly distinctive R 2 -D distributions that accurately capture the main microscopic features of the various tissues -CSF: high isotropic diffusivity D iso , low normalized diffusion anisotropy D Δ , low R 2 ; WM: low D iso , high D Δ , high R 2 ; GM: low D iso , low D Δ , high R 2 . Voxels comprising mixtures of GM, WM, and GM are in turn characterized by multimodal distributions that exhibit a linear combination of properties of the distributions from the individual components. Figure 1A displays the distribution obtained from a voxel containing both CSF and contributions from two WM tracts: the corpus callosum (CC) and the fornix. Three distinct tissue environments can be clearly discerned in the displayed distribution: an isotropic fast diffusing component attributed to CSF and two anisotropic slow diffusing components with different orientations corresponding to the WM tracts. By ascribing distribution points to one of the three bins discussed in the Methods section we were able to separate and quantify the signal contributions from distinct brain tissues. Indeed, as shown in Figure 1B, the signal fractions from the various bins follow the expected spatial distributions of WM, GM, and CSF. Figure 4 displays the ODFs computed from the components that fall within the 'thin' bin. The ODFs are displayed as directionally-coloured glyphs, superimposed on the sum of the signal fractions from the 'big' and 'thick' populations. Overall, the reconstructed ODFs are consistent with the expected WM arrangement of the healthy human brain. Major WM tracts such as the corticospinal tract (CST), the CC, and the superior longitudinal fasciculus (SLF) are easily located in the displayed figure (see arrows in Figure 4), and multiple crossings can also be discerned. The zoomed panels show that the proposed method can capture the crossings in the ventral SLF -anterior-posterior fibres with leftright fibres -and the crossings between the CST and the CC -superior-inferior fibres with left-right fibres. The dotted boxes show that three-fibre crossings present in the centrum semiovale are well captured by this technique, meaning that more than two fibre populations can be resolved. Voxels at the WM-CSF and WM-GM interfaces exhibit small-amplitude ODFs, consistent with lower signal fractions of fibrous tissue. The low amplitude of the ODF lobes found in those regions does not seem to bias their orientation; for example, CC voxels near the ventricles yield low amplitude lobes whose orientations follow the expected trend (fibres running left-right). These observations indicate that the estimated ODFs are robust to partial volume effects with CSF and that the proposed method can indeed resolve fibre orientations in heterogeneous voxels. In silico calculations show that an accurate ODF can be estimated as long as the contribution from CSF accounts for less than 75% of the total voxel-signal. Low-amplitude ODF lobes can also be found throughout cortical GM regions. These ODFs might be explained by the presence of anisotropic tissue components in cortical GM (Assaf, 2019), or interpreted as originating from low-amplitude WM partial volume effects caused by the large voxel-size used in this study. A more in-depth study than the one presented in this contribution is necessary in order to unambiguously discriminate between these two factors.

In vivo orientation-resolved R 2 -D metrics
The relaxation and diffusion features from different fibres can be investigated by using Eq. (9) to map R 2 -D metrics onto the ODF mesh and define orientation resolved means, Ê[X]. The estimated Ê [X] values are then visualized as colour-coded ODF glyphs such as the ones displayed in Figure 5, which inform on the correlations between D orientation and R 2 , D iso , or D D 2 . The displayed ODF maps capture the expected diffusion properties of healthy WM, namely a constant Ê[D iso ] ~ 1´10 -9 m 2 s -1 and a high anisotropy Ê[D D 2 ] ~ 0.7. The anisotropy metric Ê[D D 2 ] is found to be unaffected by the presence of fibre crossings (see lower right panel of Figure 5); this is in contrast to the widely popular Fractional Anisotropy (FA) metric, which is highly dependent on the degree of orientational order . Significantly lower Ê[D D 2 ] values are found at WM-GM interfaces, an observation that indicates the presence of tissue components with a lower diffusion anisotropy than that of components within pure WM voxels. Finally, we would like to note that glyphs close to ventricles do not reveal an increased D iso or decreased R 2 , thus evidencing the successful resolution of signal contributions from distinct tissue components. Focusing on the Ê[R 2 ]-coloured ODFs shown in the left side panels of Figure 5, we find a population of fibres with considerably high Ê[R 2 ] values in the midbrain region (see dashed box in the top left map of Figure 5). The fast-relaxing ODFs can be attributed to the myelinated axons that traverse the globus pallidus, an iron-rich basal ganglia structure that is characterized by particularly high R 2 values (Hasan, Walimuni, Kramer, & Narayana, 2012;Knight et al., 2015). Not accounting for their significantly different R 2 would then lead to an underestimation of the signal fraction of those high-R2 anisotropic components. Moreover, acquiring diffusion-weighted data measured at a single relatively high t E could even obscure the presence of anisotropic tissues in the pallidum.
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 6, where no clear relationship between Ê[X] and peak orientation was 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 et al., 2015;McKinnon & Jensen, 2019) 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 WM tissue. Figure 3A 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 trend observed in Figure 6 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 computed from different bootstrap ℇ b thin components. Such procedure 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). Despite the fact that no global R 2 (q) behaviour could be teased out, the proposed method allowed the detection of relaxation differences between distinct WM tracts. As shown in Figure 7, these differences are best visualized in a T 2 scale spanning a more constrained interval of values than the R 2 scale used in Figure 5. Inspection of Figure 7 reveals that both the CST and the forceps major tracts are characterized by considerably higher Ê[T 2 ] (lower Ê[R 2 ]) values. These observations are in accordance with the results of (Lampinen et al., 2020), where larger T 2 times were consistently found in the CST. The higher T 2 of the CST is also observed in voxels containing fibre crossings, with Ê[T 2 ] differences being discerned between the ODF lobes corresponding to the CST and the lobes that capture fibre populations from other tracts (see bottom right panels of Figure 7). While the exact mechanisms driving the high T 2 values found in the CST and the forceps major are still unclear, it is worth mentioning that these tracts are known to feature higher-than-average fractions of both myelin water (Coelho, Pozo, Jespersen, & Frangi, 2019) and large axons (Dell'Acqua et al., 2019). While useful for visualization purposes, the colour-coded glyphs derived in this work are however impractical for quantifying the dispersion of R 2 -D descriptors within a given ODF lobe. For example, the in silico distributions from Figure 3A demonstrate that a single fibre population may comprise a dispersion of T 2 values that cannot be accounted for by simply computing the Ê[T 2 ] value of the associated ODF peak. In this regard, we suggest using the ODFs and corresponding peaks as a guide to define additional bins in the (q,f) space and to subsequently assign the voxel-wise ℇ b thin components into the various orientation-resolved bins. Once the orientation bins have been defined and the ℇ b thin components assigned, orientation-specific statistical metrics and uncertainty measures can be estimated by exploring the variability of components within a given (q,f)-bin. An illustration of this procedure is presented in Figure 8 for an SLF voxel comprising two crossing fibres. There, the (q,f)space was divided into four quadrants centred around the extracted ODF peaks; average and dispersion measures were then calculated as the median and interquartile range of the ℇ b thin components falling within each quadrant. The procedure depicted in Figure 8 showcases the potential of using P(R 2 ,D) distributions to extract the average and variance of fibre-specific metrics. In a preliminary work , we combine the presented ODF framework with density-based clustering algorithms (Rodriguez & Laio, 2014) in order to sort ℇ b thin into different fibre populations and then calculate fibre-specific statistical metrics from the clustered P(R 2 ,D) components.

CONCLUSION
This work presents analysis protocols to estimate and visualize orientation-resolved R 2 -D metrics in the living human brain. We build on a recently developed 5D relaxation-diffusion correlation framework where sub-voxel heterogeneity is resolved with nonparametric P(R 2 ,D) distributions , and convert the retrieved distributions to ODF glyphs informing on the relaxation-diffusion features along different orientations by mapping discrete P(R 2 ,D) components to a dense mesh of (q,f) bins. Directionally-coloured ODFs estimated in such way were observed to capture fibre crossings in major WM tracts such as the CC, the CST, or the SLF. Similarly, arrays of T 2 -, R 2 -, D iso -, and D Δ 2 -coloured ODF glyphs were observed to facilitate a clean and compact visualization of the R 2 -D properties of anisotropic tissues. Maps of relaxation-coloured ODF also enabled the identification of fast-relaxing anisotropic components in the globus pallidus and the observation of long T 2 times in the CST and the forceps major.
The proposed framework relies 5D R 2 -D distributions that provide a clean 3D mapping of the signal contributions from different sub-voxel tissue environments and allow the estimation of relaxation or diffusion differences between distinct fibre populations. Moreover, the P(R 2 ,D) are retrieved from the data without the need to a priori fix signal response functions or formulating assumptions about the number of microscopic tissue components. This is in contrast with traditional (Anderson, 2005;Dell'Acqua et al., 2007;Dell'Acqua & Tournier, 2019;Tournier et al., 2007;Tournier et al., 2004) or multi-tissue  spherical deconvolution approaches, which assume a single response function for WM tissue and do not accommodate for microstructural differences across fibres. The caveat is that the proposed method hinges on signal acquisition in a high dimensional space in order to better capture the signal contrast between environments with different MR properties (Topgaard, 2019); a comprehensive sampling of such space in turn introduces acquisition times that are longer than the ones currently used in spherical deconvolution protocols. However, there is potential to reduce the scan time either by using multi-band acquisition schemes (Barth, Breuer, Koopmans, Norris, & Poser, 2016) or designing more abbreviated acquisition protocols. Recent advances in nonparametric protocol optimization (Bates, Daducci, & Caruyer, 2019;Song & Xiao, 2020) are expected to facilitate a reduction of the measured data points while keeping a good performance of the Monte Carlo inversion procedure. Protocol optimization strategies can additionally be used to maximise the angular coverage of the acquisition scheme and hopefully increase the angular resolution of the retrieved distributions (Caruyer et al., 2011). The information retrieved with the presented methodology can serve as an input for fibre tracking algorithms and used to extract individual WM pathways. If combined with tractometry frameworks (Bells et al., 2011;Chamberland et al., 2019;De Santis et al., 2014;Rheault et al., 2017;Yeatman et al., 2012), the correlations across the R2-D space would allow a comprehensive inspection of the relaxation and diffusion properties along a given WM tract. Since no universal signal response kernels are assumed, microstructural differences between tracts can be investigated and teased out. This feature is particularly promising for clinical research studies (Fornito, Zalesky, & Breakspear, 2015) where the 5D R2-D correlation framework could be used to investigate pathology induced changes along specific WM bundles.

CONFLICT OF INTERESTS
João P. de Almeida Martins, Alexis Reymbaut and Daniel Topgaard declare their status as former employee, employee, and employee/co-owner, respectively, of the private company Random Walk Imaging AB (Lund, Sweden), which holds patents related to the described method. Filip Szczepankiewicz and Daniel Topgaard are inventors on patents related to the study that are owned by Random Walk Imaging AB. The remaining authors declare no competing interests.   (6), creating an ODF glyph whose radii inform on the R 2 -D probability density along a given (q,f ) direction (middle). Following the ODF estimation, Eq. (9) is used to assign mean values of R 2 , D iso , or D Δ to each mesh point and define the colour the ODF glyph (right).

Figure 2.
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 . The displayed waveforms were inserted within a spin-echo 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 the EPI block being centred at the echo-time t E . Relaxation-weighting is enforced by varying t E while keeping constant the location of G 1 (t) and G 2 (t) relative to the 180° pulse.   Figure 1. The voxel-wise P(q,f) were computed by using Eq. (6) to map the weights of the bin-resolved discrete P(R 2 ,D) components into a 1000point spherical mesh. Here, each ODF is represented as a 3D polar plot with a local radius given by P(q,f) and colour-coded according to [R,G,B] = [µ xx ,µ yy ,µ zz ], where µ ii are the elements of the unit vector µ (q,f) (see Eq. (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 ), i.e., 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.  , 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.  thin , were then assigned to either fibre population A or fibre population B depending on their (q,f) coordinates (e.g. ℇ b 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 2 , isotropic diffusivity D iso , squared normalized diffusion anisotropy D D 2 , and orientation, thus obtaining a set of 96´5 scalars: 96 different estimates of five distinct parameters. (B) Ensemble of fibre-resolved orientations displayed on the unit sphere. The colouring of the sphere identifies the (q,f) 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 2 , D iso , and D D 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.