Reproducibility of the Standard Model of diffusion in white matter on clinical MRI systems

For Gaussian diffusion, the picture gets simplified, as the cumulant series is truncated at the level of the second-order velocity cumulant $\left<v_{i}(t)v_{j}(t^{\prime})\right>=2D_{ij}\delta(t-t^{\prime})$ defining the diffusion tensor $D_{ij}$. Hence, the signal

gets reduced to a function of the $3\times 3$ symmetric tensor $\mathbf{B}=(B_{ij})$:

where the trace of $\mathbf{B}$ is the conventional $b$-value²²2We use the “microstructure units”, i.e., $\mu$m and ms, throughout the paper. In such units, the diffusion weighting, $b=1000\,\mathrm{s/mm^{2}}=1\,\mathrm{ms/\mu m^{2}}$. In these units, free water diffusivity at body temperature $D_{\rm w}\approx 3\mathrm{\mu m^{2}/ms}$., and $\mathbf{g}(t)=\mathrm{d}{\mathbf{q}}/\mathrm{d}t$ is the diffusion gradient waveform used to generate the trajectory. Hence, for a medium (e.g., a voxel) comprised of multiple non-exchanging Gaussian-diffusion compartments (in general, anisotropic, such as in the SM), the $\mathbf{B}$-tensor fully parametrizes the measurement.

The number of nonzero eigenvalues in $\mathbf{B}$ reflects how many dimensions of the anisotropic Gaussian diffusion are being probed simultaneously. In this work, we focus on axially symmetric $\mathbf{B}$:

that are parametrized by the overall scale (the $b$-value), the unit vector $\hat{\mathbf{u}}$ along the symmetry axis, and the shape parameter $\beta$ (Fig. 2a). Compared to conventional dMRI, the extra degree of freedom $\beta$ represents the $\mathbf{B}$-tensor shape, e.g., $\beta=1$ for linear encoding (a single nonzero eigenvalue), $\beta=0$ for spherical encoding (isotropic $\mathbf{B}$-tensor), and $\beta=-\tfrac{1}{2}$ for planar encoding (two nonzero eigenvalues). Here $\delta_{ij}$ is the Kronecker symbol (the unit matrix).

The encoding time required for the gradient waveforms corresponding to a given $\mathbf{B}$-tensor depends on its shape and b-value. Figure 2b shows the minimum diffusion encoding times for $\{b,\beta\}$ combinations given specific hardware constraints (more details in Section 3.1). Isotropic weighting and large $b$-values demand more encoding time if gradient hardware constraints are kept constant.

Consider an elementary fiber segment or fiber fascicle, which is a local bundle of aligned sticks with the extra-neurite space surrounding them. The signal from such fascicle oriented along the unit vector $\hat{\mathbf{n}}$, has contributions from two axially symmetric non-exchanging Gaussian compartments aligned along $\hat{\mathbf{n}}$, Fig. 1:

• 1.

  Stick compartment, with signal fraction $f$, representing axons and possibly other elongated cells such as glial processes. Sticks are zero-radius cylinders, where diffusion occurs only along the cylinder axis with diffusivity $D_{\text{a}}$, such that the diffusion tensor is $D_{\text{a}}\,n_{i}n_{j}$.

• 2.

  Zeppelin compartment reflecting hindered diffusion in the extra-axonal space. Its diffusion tensor eigenvalues are the parallel and perpendicular diffusivities $D_{\text{e}}^{\parallel}$ and $D_{\text{e}}^{\perp}$, such that the tensor is $\Delta_{\text{e}}\,n_{i}n_{j}+D_{\text{e}}^{\perp}\delta_{ij}$, where $\Delta_{\text{e}}=D_{\text{e}}^{\parallel}-D_{\text{e}}^{\perp}$.

• 3.

  Free water compartment with signal fraction $f_{\text{w}}$ and fixed isotropic diffusivity $D_{\text{w}}=3\,\mathrm{\mu m^{2}/ms}$ is optionally added to account for CSF partial volume contributions. We will include it in our analysis.

Applying Eqs. (1) and (3) to each compartment above, and using

where $\xi=\hat{\mathbf{n}}\cdot\hat{\mathbf{u}}$ is the cosine of the angle between the symmetry axes of the kernel and of the $\mathbf{B}$-tensor, we obtain a fascicle’s response function (or a response kernel) to the measurement encoded by the $\mathbf{B}$-tensor:

Voxels contain not one but a collection of fiber segments whose orientation is given by a probability distribution on the sphere $\mathcal{P}(\hat{\mathbf{n}})$, dubbed fiber orientation distribution function (ODF). Thus, the SM signal becomes the convolution of the kernel and the ODF on the unit sphere³³3Technically, the convolution is defined on the rotation group SO(3), equivalent to the 3-dimensional unit sphere $\mathbb{S}^{3}$. However, due to the fiber fascicle’s axial symmetry, the SO(2) rotation around its axis can be factored out, and the convolution becomes over the factor group SO(3)/SO(2) equivalent to the 2-dimensional unit sphere $\mathbb{S}^{2}$, $||\hat{\mathbf{n}}||=1$ [Healy et al., 1998]. Henceforth, as in [Novikov et al., 2018], we normalize the measure on $\mathbb{S}^{2}$, $\mathrm{d}\hat{\mathbf{n}}\equiv\sin\theta\,\mathrm{d}\theta\,\mathrm{d}\phi/4\pi$, such that $\int\!\mathrm{d}\hat{\mathbf{n}}=1$.:

where $s_{0}\equiv S(\mathbf{B})|_{\mathbf{B}=0}$ is the non-weighted signal, and $\mathcal{P}(\hat{\mathbf{n}})$ is normalized to the unit probability, $\int_{\mathbb{S}^{2}}\mathcal{P}(\hat{\mathbf{n}})\,\mathrm{d}\hat{\mathbf{n}}=1$.

The SM ODF is represented via a spherical harmonic (SH) decomposition:

where $Y_{\ell m}(\hat{\mathbf{n}})$ are the SH basis functions conventionally normalized to

$p_{\ell m}$ are the SH coefficients (only even $\ell$ are nonzero due to the time-reversal symmetry of the Brownian motion and $p_{00}=\sqrt{4\pi}$), and $\ell_{\max}$ is the maximum order in the expansion (typically 4–8, depending on the SNR and the maximal $b$).

The ODF form (6) is chosen due to a number of reasons. First, the SH basis is standard for functions on a sphere, and it does not give preference to any particular functional form of the ODF (since there are currently no justifiable empirical ODF models). Second, this basis realizes the angular-radial connection in the $q$-space [Novikov et al., 2018] (the successive terms in the Taylor expansion of Eq. (5) in the powers of $q_{i_{1}}\dots q_{i_{\ell}}\sim b^{\ell/2}$ correspond to the sensitivity to SH up to order $\ell$). Third, the convolution on a sphere becomes a product in the SH basis, as discussed below. As a result, the free parameters of the SM are the compartmental diffusivities and water fractions of the kernel, and the fiber ODF coefficients $p_{\ell m}$.

When referring to the SM, unless specified otherwise, it is implied that kernel and ODF parameters are estimated directly from the data without constraints on parameters [Novikov et al., 2019]. Kernel parameters provide important tissue microstructural information, and have shown potential clinical relevance as they are sensitive to specific disease processes such as demyelination [Fieremans et al., 2012, Jelescu et al., 2016b] ($D_{\text{e}}^{\perp}$), axonal loss [Fieremans et al., 2012] ($f$) or beading [Lee et al., 2020] ($D_{\text{a}}$). Additionally, we anticipate that ODF parameters ($p_{\ell m}$) can be used for a more accurate tractography since the kernel is estimated at each voxel locally, rather than being averaged over white matter tracts as in model-free deconvolution methods [Tournier et al., 2004].

To remove the dependence on the choice of the physical basis in three-dimensional space (via $m=-\ell...\ell$), the rotational invariants of the signal and ODF are employed:

The above normalization is chosen such that: i) $p_{0}=1$, since the integral of the ODF over the sphere equals 1; ii) the remaining ODF invariants characterizing anisotropy satisfy $0\leq p_{\ell}\leq 1$. It follows that for $\ell\geq 2$ an isotropic ODF has $p_{\ell}=0$ while a delta-function on a sphere (a perfectly aligned fiber tract) has all $p_{\ell}=1$. Hence, $p_{\ell}$ are the “partial” ODF anisotropy metrics at each degree $\ell$.

From Eqs. (7) and (10), one can relate the signal rotational invariants to the kernel parameters:

This allows separating the parameter estimation in two steps $\{S(b,\beta,\hat{\mathbf{u}})\}\rightarrow\{S_{\ell}(b,\beta)\}\rightarrow\text{kernel}$, without loss of information and having to estimate only a few $p_{\ell}$ from the ODF. The above treatment generalizes earlier factorization approach [Reisert et al., 2017, Novikov et al., 2018] from LTE to arbitrary axially-symmetric $\mathbf{B}$-tensors.

A major advantage of polynomial regression over neural networks is that it is a linear optimization problem and the training can be computed much faster without the risk of local minima. Coelho et al. [2021b] provided fast analytical equations to compute not only optimal regression coefficients but also the MSE over a distribution of values for all model parameters with a given acquisition protocol and noise level. This makes MSE a good metric for comparing the performances of different protocols and thus, a better objective function than the CRB for optimal experimental design.

The metric for quantifying the goodness of a protocol was: $\text{RMSE}_{\text{obj}}=\text{RMSE}_{f}+\tfrac{1}{3}\text{RMSE}_{D_{\text{a}}}+\tfrac{1}{3}\text{RMSE}_{D_{\text{e}}^{\|}}+\tfrac{1}{3}\text{RMSE}_{D_{\text{e}}^{\perp}}+\text{RMSE}_{f_{\text{w}}}$, where each parameter is normalized by its range. This assured an even sensitivity to all parameters was kept. Individual MSE values were computed analytically as in Coelho et al. [2021b], based on the measured shells, the noise level, and the moments of the training data distribution. This accurately captures how experimental design and SNR affect our parameter estimates over a distribution of values.

We assume that the acquisition time is fixed and defined by the user. In this work we focus on 15-minute scan times. Thus, the optimization framework had to find the best way to fill the available time without violating hardware limitations such as maximum gradient amplitude and slew rate. Shells of uniformly distributed directions were assumed to be part of the optimal acquisition since it has been shown that measurements grouped into shells increase precision [Coelho et al., 2019a]. For each shell, the optimization framework selected: diffusion weighting $b$, the $\mathbf{B}$-tensor shape parameter $\beta$, number of directions, and TE.

Using the framework proposed by Sjölund et al. [2015], Szczepankiewicz et al. [2019] we generated a library of minimum encoding times as a function of the $\{b,\beta\}$ combination, see Fig. 2b. These were used to compute how much encoding time was needed for each specific $\mathbf{B}$-tensor on each scanner. The optimization handled the trade-off between adding many shorter high-SNR measurements or fewer longer low-SNR ones. The TE was constrained to be equal for all shells to factor out $T_{2}$ dependence in the analysis. Thus, the overall TE, and the SNR of the dataset, was determined by the longest diffusion waveform.

We included conventional DKI shells (LTE, $b=0-1-2\mathrm{ms/\mu m^{2}}$, [Sotiropoulos et al., 2013, Casey et al., 2018, Alfaro-Almagro et al., 2018]) as fixed into our protocols. This enables future comparisons against standard DKI-derived maps and increases the flexibility of the data analysis. The maximal $b$-value that the optimizer could explore was set to $b_{\text{max}}=10\mathrm{ms/\mu m^{2}}$ and $b_{\text{max}}=8\mathrm{ms/\mu m^{2}}$ to accommodate our two scanners described below. Stochastic optimization [Zelinka, 2004] was used to navigate the high-dimensional and non-convex protocol landscape.

On the Skyra scanner, scan(1)-rescan(2) of this protocol was acquired with TR/TE=6700/127ms (SKYRA${}_{\text{127}}^{\text{(1)}}$ and SKYRA${}_{\text{127}}^{\text{(2)}}$), while on the Prisma scanner, scan(1)-rescan(2) was acquired with TR/TE=5300/92ms (PRISMA${}_{\text{92}}^{\text{(1)}}$ and PRISMA${}_{\text{92}}^{\text{(2)}}$) in addition to one scan(3) matching the Skyra protocol with TR/TE=6700/127ms (PRISMA${}_{\text{127}}^{\text{(1)}}$). Imaging parameters: resolution: $2.0\mathrm{mm}$ isotropic, in-plane FOV: 220mm, GRAPPA and SMS acceleration factors: 2, $\text{PF}=6/8$. Subjects were taken out of scanner and repositioned between scans.

Data was subsequently processed with the DESIGNER pipeline [Ades-Aron et al., 2018]. Denoised images were corrected for Gibbs ringing artifacts [Kellner et al., 2015, Tournier et al., 2019], based on re-sampling the image using local sub-voxel shifts. These images were rigidly aligned and then corrected for eddy current distortions and subject motion simultaneously [Smith et al., 2004]. A $b=0$ image with reverse phase encoding was included for correction of EPI-induced distortions [Andersson et al., 2003].

For each kernel parameter $x$, the coefficients of variation $\mbox{COV}=\sigma/\mu$ were computed from the two repetitions $x_{1,2}$ we had for each protocol, where the mean and standard deviation were estimated as: $\hat{\mu}=\tfrac{1}{2}(x_{1}+x_{2}),\,\hat{\sigma}=\tfrac{\sqrt{\pi}}{2}|x_{1}-x_{2}|$. The, COV values were averaged over either voxels or over ROIs.

Non-LTE data contributes significantly to parameter precision. As an example, Fig. 3c shows the normalized and averaged RMSE of all model parameters after fixing the DKI shells and high-b LTE shell. Here, we move the position of the shell B in the full $\{b,\beta\}$ acquisition space. We remove shell C, and keep TE and number of directions fixed for simplicity. It can be appreciated that small modifications will not affect the protocol significantly since the RMSE landscape is very flat in their neighborhood. However, the sharp transition of the averaged RMSE as shell B moves away from LTE in Fig. 3c makes evident that after adding high-$b$ LTE data to the DKI shells, non-LTE shells are the most informative.

The COV of SM parameters are comparable to the ones from FA estimated from the DKI subset. This is encouraging since FA is known to be reproducible due to being a much simpler parameter to measure and estimate (intra-scanner voxelwise COV $\simeq 7-13\%$, ROI COV $\simeq 1-4\%$, see Fig. 6).

Interestingly, SM voxelwise COV values from noise propagation simulations agreed with the measured ones, see Fig. 6c. This agreement implies that the main factor hindering parameter reproducibility is measurement noise.

We report coefficients of variation of SM parameters around $\simeq 5-10\%$ at the voxel-level and around $\simeq 1-4\%$ for ROI means, which are actually comparable to the ones for FA (see Fig. 6). Noise propagation experiments using the estimated parameters from the WM ROIs showed our predicted COVs from simulation are similar to the ones from experiment ($\pm 50\%$), see Fig. 6c. Motion and image artifacts contribute to differences with the noise propagation predictions. Nonetheless, approximate COV estimates further support the use of the analytically predicted MSE as a quality metric for optimal protocol design.

The obtained COVs are encouraging since precise and unconstrained SM estimation has remained elusive for clinically feasible protocols. Unlike previous biophysical modeling works [Fieremans et al., 2011, Zhang et al., 2012, Jelescu et al., 2016a, Kaden et al., 2016] , we remove all parameter constraints. This makes parameter estimation more challenging but reduces biases due to arbitrary constraints and provides additional contrasts that may aid in detecting pathology. Combining optimal acquisition protocols and robust parameter estimation allowed us to obtain SM parameters from 15-minute scans, thereby enabling brain tissue microstructure mapping in clinical settings.

All SM-parameter values are found to be consistent between subjects (see Fig. 7). While ground truth values are unavailable for in vivo dMRI experiments, their values are plausible and agree well with prior studies. The axonal water fraction $f$ showed the largest values in regions with the most densely packed axons while $f_{\text{w}}$ only has large values close to the ventricles likely due to CSF partial volume. The major WM ROIs show $D_{\text{a}}$ values centered at $\sim 2.3\mathrm{\mu m^{2}/ms}$ for the TE=92ms data, in agreement with measurements using more extensive acquisitions involving high diffusion weighted PTE done by Dhital et al. [2019], and with high b LTE [Veraart et al., 2019a], yet slightly lower than those reported by [Howard et al., 2020, Nilsson et al., 2021]. Furthermore, $D_{\text{a}}>D_{\text{e}}^{||}$ in the majority of voxels, agreeing with gadolinium-based contrast experiments in the rat corpus callosum [Kunz et al., 2018], though this observation depends on the ROI and both values are generally close to each other.

Due to duty cycle gradient heating, LTE high-b shells were acquired at $b=5.5\,\mathrm{ms/\mu m^{2}}$ rather than $b=8\,\mathrm{ms/\mu m^{2}}$ for the Prisma and at $b=5\,\mathrm{ms/\mu m^{2}}$ rather than $b=7\,\mathrm{ms/\mu m^{2}}$ for the Skyra, as initially suggested by the protocol optimization to avoid increasing the TE (see Fig. 3a). Noise propagation experiments showed such change did not significantly deteriorate the quality of the protocol. This happens because at high-b the location of shell A is already independent from the rest and the MSE becomes very flat around it, see a similar situation for shell B in Fig. 3c.

Overall, the protocol optimization is a high-dimensional problem which has multiple local minima. Our approach showed robust estimations of global optima in a toy function of similar dimensionality and complexity. We do not explore $\mathbf{B}$-tensor shapes without axial symmetry (inner part of the triangle in Fig. 2a). Releasing this constraint would increase the space of acquisitions but not necessarily increase precision in the parameter estimation, as shown previously for the cumulant expansion parameters Coelho et al. [2019b]. Therefore, we constrained $\mathbf{B}$-tensors to be axially symmetric, reducing the dimensionality of the optimization problem. Future work will explore reproducibility of optimal acquisitions with varying TE to simultaneously capture compartmental $T_{2}$ values, as it was proposed by Veraart et al. [2018], McKinnon and Jensen [2019] and later assessed by Lampinen et al. [2020].