Synopsis

Keywords: Microstructure, Diffusion/other diffusion imaging techniques

Microstructure imaging with diffusion MRI relies on non-linear fitting of a biophysical tissue model to the data, often at low signal-to-noise ratio. In this work, we derive a new rotation-invariant feature set for microstructure mapping based on a rank-1 decomposition of the multi-shell diffusion MRI signal in spherical harmonics. Simulations show that using this feature set avoids non-central-χ bias that is present at low SNR in parameter estimation based on conventional rotation-invariants. Results in human brain imaging data acquired with free waveform diffusion encoding show robust parameter estimates across white matter.

Introduction

Diffusion MRI (dMRI) is sensitive to tissue microstructure on a scale 2 – 3 orders of magnitude below the image resolution. Microstructure models of the dMRI signal response aim to identify more specific parameters of the tissue composition and how it is affected in pathology ¹. However, robust microstructure imaging critically relies on accurate parameter estimators.

In brain white matter, the dMRI signal is well represented as the spherical convolution of a single-fibre response function $$$H(\theta | B, \mathbf{x})$$$ and an orientation distribution function (ODF) $$$P(\mathbf{u})$$$ ^2-4. The response function is often modelled as an intra-axonal “stick” compartment with fraction $$$f$$$ and axial diffusivity $$$D_i$$$ and an extra-cellular compartment with axial and radial diffusivities $$$D_e$$$ and $$$D_r$$$ ^5,6. Microstructure mapping aims to estimate the vector of microstructure parameters $$$\mathbf{x}$$$ in each voxel for which this model best fits the acquired data.

Methods

Theory. In the basis of spherical harmonics (SH), spherical convolution of the microstructure kernel and ODF can be expressed as the outer product of a vector $$$\mathbf{h}_\ell = [ \ldots h_{\ell,B} \ldots]$$$ (indexed in $$$B$$$, the free-waveform diffusion encoding) and $$$\mathbf{p}_\ell = [ \ldots p_\ell^m \ldots]$$$ (indexed in $$$m$$$) in each harmonic band $$$\ell$$$ ⁴. Using this notation and assuming Gaussian noise, microstructure mapping solves a non-linear least-squares problem:
$$ \min_{\mathbf{x}} \sum_\ell \| S_\ell - \mathbf{h}_\ell(\mathbf{x}) \mathbf{p}_\ell^\top \|^2 \quad, $$
where $$$S_\ell$$$ are matrices that structure the SH representation of dMRI signal per shell $$$B$$$ along the rows. Since, without noise, these matrices are of rank=1, we can approximate $$$\mathbf{p}_\ell = p_\ell \mathbf{v}_\ell$$$ where $$$\mathbf{v}_\ell$$$ is the leading right singular vector of $$$S_\ell \approx \mathbf{u}_\ell \sigma_\ell \mathbf{v}_\ell^\top$$$. Substituting into Eq. (1) leads to the equivalent cost function:
$$ \min_{\mathbf{x}} \sum_\ell \| \mathbf{r}_\ell - p_\ell \mathbf{h}_\ell(\mathbf{x}) \|^2 \quad, $$
where $$$\mathbf{r}_\ell = \sigma_\ell \mathbf{u}_\ell$$$ are voxel-level rank-1 SHARD features ⁴.

This equation corresponds to rotation-invariant microstructure fitting ⁶ except for the feature set that is used. Prior work used the SH power spectrum (also known as RISH-features), i.e., the norm of each row of $$$S_\ell$$$. These features are affected by non-central-χ bias at low SNR, whereas the derivation of Eq. (2) holds at any SNR.

Data. Data were collected in a healthy volunteer on a Philips 3T Achieva MRI using a 32-channel head coil. TE = 90ms; TR = 6000ms; slice thickness = 2.5mm; in-plane resolution = 2.5mm; MB = 2; SENSE = 2. The diffusion encoding comprises 3 shells with linear encoding at b=1000, 2000 and 5000 s/mm², and 2 shells with b_Δ = 0.6 at b=1000 and 2000 s/mm², following ⁷. Each shell consists of 60 gradient directions, in addition to 16 b=0 images. Data were preprocessed with denoising ⁸, field unwarping ⁹, and motion correction ¹⁰.

Parameter fitting. We implemented rotation-invariant model fitting ⁶, using constrained Levenberg-Marquardt (LM) least-squares optimization with variable projection for the linear parameters $$$p_\ell$$$ ¹¹. The Jacobian is calculated analytically. Iterative optimization is initialised using a Baysian Monte-Carlo estimator. All experiments are run at $$$\ell_{max} = 8$$$. The run time for full brain analysis is 1min 4sec on an Apple M1 Pro CPU.

Results

Simulations. Figure 1 plots the parameter landscape for a “true” value (black cross) set as the mean in white matter reported in ⁷ and using the chosen diffusion encoding scheme. We then simulated 100 realisations of Gaussian noise at SNR=30 and, initialised at the ground truth, and ran LM fitting using either RISH (orange) or SHARD (green) features. While the latter cluster around the ground truth, the RISH-based optimization shows significant bias in all parameters (univariate t-test to ground truth).

Figure 2 repeats this experiment for a broad range of simulated ground truth values. While the RISH-based estimator introduces bias and even splits into two branches, the SHARD-based estimator is centred around a line of unit slope.

In vivo data. Figure 3 shows estimated parameter maps in human brain data. We find that usually $$$D_e > D_i$$$. Differences between RISH-based and SHARD-based estimators are less outspoken than the simulations predicted but arise in the extra-axonal diffusivity in the superior longitudinal fasciculus and the pallidum (Table 1).

Discussion and Conclusion

We presented a rederivation of rotation-invariant microstructure mapping with a subtle conceptual difference: instead of defining rotation-invariant features from the outset as a means to factor out the ODF, we showed that substituting the best model-free guess of the ODF results in a similarly compact feature set.

The simulation results highlight a potential source of bias in RISH-based microstructure mapping. This non-central-χ bias is dominant at low SNR, i.e., at high b-values and at high SH order $$$\ell$$$. The simulations also illustrate that the feature set derived in Eq. (2) is not affected by this bias. Results in human brain data are more subtle, possibly due to non-Gaussian noise effects.

Despite clear improvements in parameter bias, the large variance of the fitted parameters (especially $$$D_e$$$ and $$$D_i$$$) in dMRI data at realistic SNR remains a considerable challenge for maximum-likelihood estimation. Future work can design Bayesian estimators ⁵ based on the proposed feature set.

Acknowledgements

DC is supported by the Flemish Research Foundation (FWO; grant number 12ZV420N).