Toward nonparametric diffusion‐T1 characterization of crossing fibers in the human brain

Purpose To estimate T1 for each distinct fiber population within voxels containing multiple brain tissue types. Methods A diffusion‐T1 correlation experiment was carried out in an in vivo human brain using tensor‐valued diffusion encoding and multiple repetition times. The acquired data were inverted using a Monte Carlo algorithm that retrieves nonparametric distributions P(D,R1) of diffusion tensors and longitudinal relaxation rates R1=1/T1. Orientation distribution functions (ODFs) of the highly anisotropic components of P(D,R1) were defined to visualize orientation‐specific diffusion‐relaxation properties. Finally, Monte Carlo density‐peak clustering (MC‐DPC) was performed to quantify fiber‐specific features and investigate microstructural differences between white matter fiber bundles. Results Parameter maps corresponding to P(D,R1)’s statistical descriptors were obtained, exhibiting the expected R1 contrast between brain tissue types. Our ODFs recovered local orientations consistent with the known anatomy and indicated differences in R1 between major crossing fiber bundles. These differences, confirmed by MC‐DPC, were in qualitative agreement with previous model‐based works but seem biased by the limitations of our current experimental setup. Conclusions Our Monte Carlo framework enables the nonparametric estimation of fiber‐specific diffusion‐T1 features, thereby showing potential for characterizing developmental or pathological changes in T1 within a given fiber bundle, and for investigating interbundle T1 differences.


S1 ORIENTATION DISTRIBUTION FUNCTIONS
For each bootstrap solution n b (with 1 ≤ n b ≤ N b = 96), we considered the voxel-wise discrete ensemble of components belonging to the thin bin, and computed an ODF Pn b (θ mesh , φ mesh ) on a dense spherical mesh {(θ mesh , φ mesh )} by convolving the discrete set of components of Equation S1 with a Watson distribution kernel 4,5 as where u i ≡ (θ i , φ i ) is the unit vector giving the orientation of component i, µ(θ mesh , φ mesh ) ≡ (θ mesh , φ mesh ) is the unit vector associated with a spherical-mesh point, "·" denotes the scalar product, and κ is the concentration parameter that regulates the amount of orientation dispersion around u i in the Watson kernel. A final voxel-wise ODF P (θ mesh , φ mesh ) was calculated as the median of the per-bootstrap ODFs: The purpose of the Watson kernel is to smoothly map the discrete set of components onto the nearest nodes of the spherical mesh {(θ mesh , φ mesh )} without inducing any significant peak broadening larger than the distance between two nearest-neighboring mesh nodes. In this work, we considered a 1000-point uniform spherical mesh and set the concentration parameter of the Watson kernel to κ = 14.9, corresponding in the small-angle limit to a "spherical Gaussian" with standard deviation equal to 1 2 and R 1,i values enable the computation of orientation-specific diffusion-relaxation means: 1,2 with χ ≡ D iso , D 2 ∆ , R 1 . After isolating the peak orientations {(θ peak , φ peak )} where the median ODF P (θ mesh , φ mesh ) is locally maximized, ODF-peak measures of χ ≡ D iso , D 2 ∆ , R 1 can be obtained by computing Equation S4 at any given peak orientation.

MONTE-CARLO DENSITY-PEAK CLUSTERING
Let us review the steps through which MC-DPC estimates fiber-specific properties from the output of the Monte-Carlo inversion algorithm. First, ) and delineates Nc clusters in its orientation subspace using DPC with data-point density and outlier detection altered to account for the weights w i of the retrieved thin-bin components. 3 An initial number of clusters Nc is automatically set by the number of voxel-wise ODF peaks output by our Monte-Carlo framework, 1,2 but may be reduced during the MC-DPC procedure following a filtering approach detailed in Ref. [3]. Second, MC-DPC computes orientation-resolved statistics across bootstrap solutions by separately classifying each per-bootstrap ensemble of thin-bin solutions E thin n b into Nc ensembles E thin n b ,nc (with 1 ≤ nc ≤ Nc), each containing the thin-bin solutions of bootstrap solution n b that belong to an estimated cluster nc. It then averages the properties of the solutions within each ensemble E thin n b ,nc independently, yielding the orientation-resolved means where (x , y, z ) are the Cartesian coordinates associated with a component orientation (θ, φ). The total weight associated with each orientation-resolved mean is given bẙ The short-hand notations "E[χ]" and "ẘ " are now used for simplicity to describe the collection of orientation-resolved meansE[χ]n b ,nc and weights wn b ,nc originating from all bootstrap solutions n b and all clusters nc. Finally, one can extract the median and interquartile range of the orientationresolved meansE[χ] and weightsẘ across bootstrap solutions.
From a pure orientation standpoint, an equivalent of a median cluster orientation can be computed as the geometric median orientation u Med,nc ≡ (θ Med,nc , φ Med,nc ) of each cluster-specific collection of mean orientations, {E[u]n b ,nc } 1≤n b ≤N b , via the argument minimization of the following function dependent on arbitrary unit orientations v: using the antipodally symmetric angular distance

S3 NUMERICAL DATA
An in silico evaluation of MC-DPC was performed by designing a set of components simplicity. Orders of magnitude for D iso , D ∆ and T 1 were inspired by Ref. [6], giving the following compartments: • one isotropic component, D iso = 2 µm 2 /ms, T 1 = 4 s, w iso = 0.1.
• one fiber along z (θ = 0, φ = 0), total weight = (1 − w iso ) × 0.41, The ground-truth signals associated with these systems were computed using the same inversion kernel and acquisition scheme as those detailed in the main body of this paper. This calculation is in agreement with the conventional procedure for testing 1D, 2D or 4D Laplace inversion algorithms, where the ground-truth signal is calculated with the same kernel as the inversion. [7][8][9][10][11] Rician noise was then added to the ground-truth signals according to where Sgt is a ground-truth signal, S is the corresponding noisy signal, ν and ν denote random numbers drawn from a normal distribution with zero mean and unit standard deviation, and SNR is the desired signal-to-noise ratio. Three SNRs were considered: the SNR = 40 of our in vivo dataset, and higher SNR values, namely SNR = 70 and SNR = 90. For each of the aforementioned SNRs, the signals obtained from Equation S9 were inverted using the 5D Monte-Carlo inversion detailed in the main body of the paper with either N b = 50, N b = 75 or N b = 100 bootstrap solutions.
Indeed, varying N b may be relevant when validating MC-DPC because the number of data points that serve as its input is roughly proportional to N b . 3 Finally, MC-DPC was run on the solutions retrieved from each of these inversions and noise levels to extract orientation-resolved means and weights (see Equations S5 and S6). Figure S1 displays the sub-voxel orientations retrieved for the in silico data described in Section S3 using both the ODFs of Section S1 and MC-DPC from Section S2, and quantifies the angular deviation∆β of each cluster geometric median orientation (see Equation S7) compared to the closest ground-truth fiber orientation. Figure S2 quantifies the orientation-resolved meansE[χ] (see Equation S5) for the same in silico data. ODF-peak information is also displayed on Figures S1 and S2 for comparison.

S4 NUMERICAL EVALUATION
Within Figure S1, the orientational spread of the MC-DPC clusters informs on the precision of the underlying Monte-Carlo signal inversion algorithm. Clusters become more orientationally dispersed as SNR decreases at constant N b , which is associated with the loss of precision of the 5D Monte-Carlo inversion with reduced SNR. At constant SNR, the orientational clusters retrieved from MC-DPC are rather unaffected upon changing N b . As for the angular deviation∆β, informing on the accuracy of the underlying Monte-Carlo signal inversion algorithm, it increases with decreasing SNR, especially at SNR = 40. Compared to ODF-peak orientations, cluster geometric median orientations seem to be more accurate as they typically yield smaller values of∆β, probably due to the fact that MC-DPC clusters are not bound to a discrete spherical mesh. However, ODF-peak orientations appear to be more resilient to noise.
Regarding Figure S2, the estimations of the orientation-resolved means and weights tend to suffer from the same loss of accuracy and precision mentioned for the orientational information in Figure S1 as SNR decreases at constant N b . Estimation performances are not affected by However, the SNR = 40 case, i.e. that closest to the in vivo study presented in the main body of the paper, seems to poorly estimateE[D iso ] and, importantly,E[T 1 ] for certain fiber populations. Similar discrepancies are reported forE[T 1 ] in the in vivo study presented in the main body of the paper. As for ODF-peak metrics, they generally agree with the median orientation-resolved means at all SNR levels.

FIGURE S1
Sub-voxel orientations retrieved for the in silico data described in Section S3 using the Monte-Carlo inversion for various numbers N b of bootstrap solutions and various SNR levels. While the ODFs were obtained via the process detailed in Section S1, the orientational clusters, here represented on the unit sphere, were extracted via MC-DPC according to Section S2.∆β denotes the angular deviation, computed for a given orientational cluster as the shortest angle between either the cluster geometric median orientation (circles, see Equation S7) or the corresponding ODF peak (squares, see Section S1), and the closest ground-truth anisotropic component orientation. The color mapped onto the ODF codes for local orientation according to [red, green, blue] ≡ [|x |, |y|, |z |]/max([|x |, |y|, |z |]). As for the clusters, while opacity codes for the weight of the intracluster averaged components (see Equation S6), color codes for the geometric median orientation of each cluster (see Equation S7). The conditions of the in vivo study presented in the main body of the paper are closest to the case (N b = 100, SNR = 40).
While ground-truth is shown as horizontal lines, the circles and whiskers represent the medians and interquartile ranges of the orientation-resolved means across bootstrap solutions, respectively. Squares correspond to the estimated ODF-peak metrics. Colors match those of the orientational clusters/ODF peaks presented in Figure S1. In the rightmost panels, cluster weightsẘ were normalized so that the sum of all median weights across clusters equals one. Their ODF-peak equivalents were simply obtained by taking the mesh-projected component weights (i.e. ODF radii) along the peaks of a given ODF (see Section S1). These ODF-peak weights were then normalized to sum up to one, for easier comparison with normalized cluster weights. The conditions of the in vivo study presented in the main body of the paper are closest to the case (N b = 100, SNR = 40).