Introduction

Quantitative MRI (qMRI) allows the characterisation of the biological structures in the human brain⁶. For example, apparent transverse R₂*-weighted gradient-echo (GRE) MR images have been shown to be sensitive to both brain microstructure (e.g. fibres's myelination⁴) and the mean angular orientation (θ_μ) of fibre pathways relative to the main magnetic field, B₀^8,9. We recently introduced an inverse model that allows separation of the orientation-independent part of R₂* and an orientation-dependent higher-order term from a single multi-echo GRE measurement⁵. That model was based on the biophysical hollow cylinder fibre model (HCFM)⁸ assuming fully parallel fibres. Here we propose an extension of our model⁵ which heuristically includes the effect of fibre dispersion in the R₂* estimation. To investigate the effects of dispersion on the GRE signal and the performance of the signal models, simulations were performed by combining the ensemble-averaged signal (EAS) model from the HCFM with angular fibre dispersion modelled by the Watson distribution².

Methods

Forward model: To study the effect of dispersion on the (R₂*) signal decay, an EAS with dispersedly oriented cylinders (defined by κ, Fig.1B) around θ_μ with respect to B₀ was simulated. In this experiment, we considered the idealized scenario⁷ where the myelin water compartment was neglected (S_M≈0, R₂*≈0), i.e. the simulated signal came from 1500 cylinders with intra-(S_A) and extra-(S_E) axonal compartments (Fig.1A). The averaged complex signal was sampled 5000 times with additive Gaussian complex noise to achieve a signal-to-noise ratio (SNR) of 100 (similar experimental SNR⁵) and infinity (ground truth)³. Simulation sampling scheme: θ_μ:[1:1:90]°, κ:[0.001:1.0:20,30:10:100,1e8 (perfectly parallel)] and echo time (TE):[0:0.1:40]ms.
Inverse model and experimental conditions: Ten time-points evenly selected between 3-36 ms from simulated data were chosen to replicate experimental conditions⁵ and fulfil the static dephasing regime⁸. These signals were then fitted to 4 different models $$$m={M1,M2,M3,M4}$$$ of the form:
$$ln(|S_{WB}|)\approx\beta_0^m+\beta_1^mTE+\alpha_2^m{TE}^2$$
where $$$\alpha_2^m$$$ is defined as follows: $$$\alpha_2^{M1}=0$$$ (classical model, model 1), $$$\alpha_2^{M2}=\beta_2^{M2}$$$ is dependent on $$$\theta$$$ (model 2), $$$\alpha_2^{M3}(\theta)=\beta_2^{M3}\sin(\theta)^4$$$ (model 3) and $$$\alpha_2^{M4}(\theta)=\beta_2^{M4} <sin(\theta)^4>_{\theta^\prime,\kappa}$$$ (model 4). For the latter expression, the averaging of the angular component was performed using the Watson distribution, giving:

$$<\sin(\theta)^4>_{\theta',\kappa}=<1-2\cos(\theta)^2+\cos(\theta)^4>_{\theta',\kappa}$$
$$=(1-2\tau+\gamma)-(1+6\tau-5\gamma)\sin(\theta_\mu)^2+(3/8+35\gamma/8-15\tau/4)\sin(\theta_\mu)^4$$
$$\tau(\kappa)=\frac{2\kappa\exp(\kappa)-\sqrt{\pi\kappa}Erfi(\sqrt{\kappa})}{2\sqrt{\pi\kappa^3}Erfi(\sqrt{\kappa})},\gamma(\kappa)=\frac{3\sqrt{\pi}Erfi(\sqrt{\kappa}) + 2\sqrt{\kappa}\exp(\kappa)(2\kappa-3)}{4\sqrt{\pi}\kappa^2Erfi(\sqrt{\kappa})};Erfi(\sqrt{\kappa})=\frac{2}{\sqrt{\kappa}}\int_0^\sqrt{\kappa}exp(t^2)dt$$

where κ is the dispersion parameter and θ_μ is the angle between $$$\vec{B_0}$$$ and the mean orientation of the Watson distribution $$$\vec{\mu}$$$ (Fig.1B). The fitting of the model parameters, by given κ-θ_μ, was performed by linear regression and stabilised with a zeroth-order Tikhonov regularization for model 3-4 (of the form of λcos(θ)⁴ by targetting only the second-order component of the signal model). To estimate the optimal regularization parameter λ₀ for SNR=100, we minimized the difference between model parameters at SNR=∞ and estimated parameters for a range of λ=1e-4:1e-4:1e-1 (λ₀ of 0.01,0.02 for model 3-4, respectively). To compare the effect of SNR in the estimation of the model parameters for each model at different κ-θ_μ, the coefficient of variation (CV) of the estimated β’s across 5000 noise configurations were calculated.

Results and discussion

The simulated GRE signals in the case of dispersion are presented in Fig.2. The signal decay lost its angular dependency for large dispersion factors (solid lines for κ≈0) whereas it resembled the ideal case of parallel cylinders (dashed lines) for low dispersion factors (solid lines for κ>20). The parameters of model 1-2 behaved as previously reported⁵, demonstrating that the orientation dependence of β₁^M1 is transferred to β₂^M2 (Fig.3A,B). β₂^M3 values, at higher angles, were similar to β₂^M2, but for smaller angles (between 15°-30°) the differences increased exponentially; this divergence was down-weighted by the regularization leading to convergence towards zero (Fig.3C). For β₂^M4, the exponential increase for low angles for β₂^M3 was compensated, showing stable values for all the angles >20° (Fig.3D). The resulting estimation of the β parameters for all the models for the experimental SNR=100 is shown in Fig.4. β₁^{M2, M3 and M4} showed a 5x increase in comparison to β₁^M1. In addition, β₂^{M3 and M4} became greatly sensitive to noise for angles <45° across all κ (CV>1, see Table 1). Nevertheless, there was no bias introduced in the mean estimation of all the parameters in all the models, i.e. the estimated β’s are close to the ground truth (SNR=∞, Fig.3).

Conclusion

In this work, we introduced a forward model to simulate the effect of fibre dispersion on the angle-dependent GRE signal decay in different SNR conditions. Moreover, we proposed a new inverse model that accounts for this effect on the estimated parameters and compared it with published models⁵ for an experimentally-relevant SNR regime. We found noise enhancement in the estimated parameters of the second-order models (M2,M3,M4) as compared to M1, which can be explained by the increased condition number in the inverse model¹ when introducing the second-order decay parameter β₂. Interestingly, we found for β₂ noise enhancement for highly aligned fibres for angles <45°, which might be explained by the vanishing quadratic signal-contribution for low angles. Nevertheless, models 1-2 proved to be still valid even for moderate dispersion. Future work will assess the errors made by neglecting the myelin compartment and the position dependence of the frequency offsets, as already simulated¹⁰.

Acknowledgements

This work was supported by the German Research Foundation (DFG Priority Program 2041 "Computational Connectomics”, [AL 1156/2-1;GE 2967/1-1; MO 2397/5-1; MO 2249/3–1], by the Emmy Noether Stipend: MO 2397/4-1) and by the BMBF (01EW1711A and B) in the framework of ERA-NET NEURON and the Forschungszentrums Medizintechnik Hamburg (fmthh; grant 01fmthh2017). The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013) / ERC grant agreement n° 616905.