Investigating the Benefits of Incorporating Higher Order Spherical Harmonics in Axon Diameter Measurements
Qiuyun Fan1,2, Aapo A Nummenmaa1,2, Qiyuan Tian1,2, Ned A Ohringer1,2, Thomas Witzel1,2, Lawrence L Wald1,2,3, Bruce R Rosen1,2,3, and Susie Y Huang1,2,3

1Radiology, Massachusetts General Hospital, Charlestown, MA, United States, 2Harvard Medical School, Boston, MA, United States, 3Harvard-MIT Division of Health Sciences and Technology, Massachusetts Institute of Technology, Cambridge, MA, United States


Separating out the scalar and orientation-dependent components of the diffusion MRI signal offers the possibility of increasing sensitivity to microscopic tissue features unconfounded by the fiber orientation. Recent approaches to estimating apparent axon diameter in white matter have employed spherical averaging to avoid the confounding effects of fiber crossings and dispersion at the expense of losing sensitivity to effective compartment size. Here, we investigate the feasibility and benefits of incorporating higher-order spherical harmonic (SH) components into a rotationally invariant axon diameter estimation framework and demonstrate improved precision of axon diameter estimation in the in vivo human brain.


Approaches for estimating microstructural tissue properties by diffusion MRI have evolved over the past decade, with increasing recognition that the scalar parameters describing biophysical features of tissue, e.g., volume fraction, diffusivity, and effective compartment size, and the orientation-dependent parameters describing the fiber orientation distribution function, carry separate and complementary information 1-6. Recent approaches to estimating apparent axon diameter in white matter using high b-value data have employed spherical averaging to avoid the confounding effects of fiber crossings and orientation dispersion 7,16. While spherical averaging enables estimation of biophysical parameters independent of crossing fibers, averaging the diffusion MRI signal over all directions reduces sensitivity to apparent axon diameter by reducing the contribution of the signal perpendicular to the principal fiber direction. The goal of this study is to investigate the feasibility and benefits of incorporating higher-order spherical harmonic (SH) components into the axon diameter estimation framework.


White matter is modeled as composed of two compartments: restricted diffusion in cylindrical axons modeled using the Gaussian phase distribution approximation8, and extra-axonal hindered diffusion modeled as Gaussian diffusion with parallel diffusivity of 1.7x10-3mm2/s. The signal was projected to the SH space (Figure 1), so that $$$s_{lm}=c_{l}\cdot p_{lm}$$$, where $$$s_{lm}$$$, $$$c_{l}$$$, and $$$p_{lm}$$$ are the SH coefficients for the diffusion signal, response kernel (Figure 2) and orientation distribution function (ODF). Then both sides of the equation were divided by a normalization factor $$$\sum_{j=1}^Ts_{lm}=\sum_{j=1}^Tc_{l}\cdot p_{lm}$$$ over all $$$T$$$ experiments, so that the $$$p_{lm}$$$ on the right side is canceled out, yielding the relationship between the normalized signal and the normalized response kernel, i.e., $$$ns_{l}$$$ = $$$nc_{l}\equiv \frac{c_{l}}{\sum_{j=1}^Tc_{l}}$$$, an ODF independent form response kernel. Specifically, $$$c_{l}=f_{r}\cdot c_{r}(a)+f_{h}\cdot c_{h}(D_{h})$$$, where $$$f_{r}+f_{h}=1$$$, $$$c_{r}$$$ is the SH representation of the van Galderen model8,9, and $$$c_{h}$$$ is the SH representation of the tensor model10.


A healthy subject was scanned on the 3T Connectome scanner with 300mT/m maximum gradient strength using a custom-made 64-channel head coil11. Real-valued diffusion data was acquired to avoid buildup of the noise floor12. Sagittal 2-mm isotropic resolution diffusion-weighted spin-echo EPI images were acquired with whole brain coverage. The following parameters were used: TR/TE=4000/77ms, δ=8ms, Δ=19/49ms, 8 diffusion gradient strengths linearly spaced from 30-290mT/m per Δ, 32-64 diffusion directions, parallel imaging (R=2) and simultaneous multislice (MB=2). Diffusion data were corrected for susceptibility and eddy current distortions using the TOPUP13 and EDDY14,15 tool in FSL.

The normalized kernel $$$c_{l}$$$ was evaluated in the 3D parameter space X=(a, Dh, fr), and the best fit to the model was found by searching on the grid by minimizing the “energy” function1 (Figure 3), i.e., $$$\widetilde{x}=\arg_{x \in X} min\sum_j^T\sum_{l=-L}^L (ns_{l}-nc_{l} (x))^{2}$$$

Simulation data was generated by adding 100 samples of noise at SNR=20. Voxel-wise fitting for axon diameter a, restricted fraction fr, and hindered diffusivity Dh was performed.


Simulation data shows that by combining L=0 and L=2 signal components, the standard deviation in axon diameter estimates decrease compared to using L=0 components alone (Figure 4). In the in vivo data, the feature of larger apparent axon diameter is more pronounced in the corticospinal tract in the L0+L2 representation by incorporating higher-order SH components (Figure 5).


We propose a method for measurement of compartment size and volume fraction by supplementing spherical averaging with higher-order SH components. While spherical averaging is robust to fiber crossings and image noise, only 0th order signal components were used with higher frequency signal components being wasted. Similar approaches have been pursued previously to make use of the additional information available in multi-shell data 1,3. Our results indicate that it is feasible to incorporate the higher-order SH in measuring compartment size by transforming the signal equation into an ODF invariant form. This approach can be applied to whole brain analyses, with the potential to improve the precision of axon diameter estimates.


This work was funded by an NIH Blueprint for Neuroscience Research Grant: U01MH093765, as well as NIH funding from NCRR P41EB015896, NIBIB R01EB006847, NIBIB R00EB015445, NINDS K23NS096056, NINDS K23NS078044, NIH/NCRR/NIBIB P41EB015896 and Instrumentation Grants S10-RR023401, S10-RR023043, and S10-RR019307. Funding support was also received from the National Multiple Sclerosis Society, the American Heart Association Postdoctoral Fellowship Award (17POST33670452), a Radiological Sciences of North America Research Resident Grant, the Conrad N. Hilton Foundation and the MGH Executive Committee on Research Fund for Medical Discovery Fellowship Award.


1. Novikov, D.S., Veraart, J., Jelescu, I.O. & Fieremans, E. Rotationally-invariant mapping of scalar and orientational metrics of neuronal microstructure with diffusion MRI. NeuroImage 174, 518-538 (2018).

2. Veraart, J., Fieremans, E. & Novikov, D. Quantifying neuronal microstructure integrity with TE dependent Diffusion Imaging (TEdDI) in Proc. ISMRM, Vol. 0836 (Honolulu, HI, USA, 2017).

3. Reisert, M., Kellner, E., Dhital, B., Hennig, J. & Kiselev, V.G. Disentangling micro from mesostructure by diffusion MRI: A Bayesian approach. NeuroImage 147, 964-975 (2017).

4. Kaden, E., Kelm, N.D., Carson, R.P., Does, M.D. & Alexander, D.C. Multi-compartment microscopic diffusion imaging. Neuroimage 139, 346-359 (2016).

5. Kaden, E., Kruggel, F. & Alexander, D.C. Quantitative mapping of the per-axon diffusion coefficients in brain white matter. Magn Reson Med 75, 1752-1763 (2016).

6. Jespersen, S.N., et al. Neurite density from magnetic resonance diffusion measurements at ultrahigh field: Comparison with light microscopy and electron microscopy. NeuroImage 49, 205-216 (2010).

7. Fan, Q., et al. Axon Diameter Mapping Independent of Crossing Structures using Spherical Mean Technique. in Proc. ISMRM 5244 (Paris, France, 2018).

8. van Gelderen, P., Moonen, C.T. & Duyn, J.H. Susceptibility insensitive single shot MRI combining BURST and multiple spin echoes. Magn Reson Med 33, 439-442 (1995).

9. Panagiotaki, E., et al. Compartment models of the diffusion MR signal in brain white matter: a taxonomy and comparison. Neuroimage 59, 2241-2254 (2012).

10. Anderson, A.W. Measurement of fiber orientation distributions using high angular resolution diffusion imaging. Magn Reson Med 54, 1194-1206 (2005).

11. Setsompop, K., et al. Pushing the limits of in vivo diffusion MRI for the Human Connectome Project. NeuroImage 80, 220-233 (2013).

12. Eichner, C., et al. Real diffusion-weighted MRI enabling true signal averaging and increased diffusion contrast. Neuroimage 122, 373-384 (2015).

13. Andersson, J.L., Skare, S. & Ashburner, J. How to correct susceptibility distortions in spin-echo echo-planar images: application to diffusion tensor imaging. Neuroimage 20, 870-888 (2003).

14. Andersson, J.L. & Sotiropoulos, S.N. An integrated approach to correction for off-resonance effects and subject movement in diffusion MR imaging. NeuroImage 125, 1063-1078 (2016).

15. Andersson, J.L.R., Graham, M.S., Zsoldos, E. & Sotiropoulos, S.N. Incorporating outlier detection and replacement into a non-parametric framework for movement and distortion correction of diffusion MR images. NeuroImage 141, 556-572 (2016).

16. Veraart J, Fieremans E, Novikov DS. On the scaling behavior of water diffusion in human brain white matter. NeuroImage. 2018 Oct 4;185:379-87.


Figure 1. Illustration of the signal representation using spherical harmonics. The white matter signal (S) can be modeled as a mixture of restricted (Sr) and hindered (Sh) diffusion signal, which can be represented by the spherical harmonics (Ylm) and corresponding coefficients (slm). When the fiber axis is aligned with the primary axis of the coordinate system, slm=0 for all m≠0, so that the signal can be fully represented by sl0, l=0,2,4,… s00 and s20 as a function of gradient strength are plotted for the restricted and hindered compartments. Note that the L=0 component is identical to the spherical mean signal.

Figure 2. The response kernel in the spherical harmonic space. The restricted diffusion response kernel (cr) as a function of diameter is shown on the left (a-d), and the hindered response kernel (ch) as a function of hindered diffusivity is shown on the right (e-h).

Figure 3. Exemplary “energy” function for different spherical harmonic orders L. The “energy” function was evaluated for the ground truth parameter {a, Dh}={6μm,10-9 m2/s} (marked with the white cross) and its neighborhood in the parameter space X = {a, Dh}. The global minimum of the “energy” appears at the ground truth value for both L=0 and L=2 components, each of which shows a slightly different pattern. By combining the two components (L0+L2), the “low-energy” region shrinks, indicating that incorporating complementary contrasts originating from different orders is helpful for finding the global minimum.

Figure 4. Simulation results of estimated axon diameters and restricted volume fractions. The in vivo imaging protocol was used to generate the noise-free data (keeping fr=0.7 for the left figure and keeping diameter=6 μm for the right), and 100 samples of noise were added with an SNR=20. The mean value for the 100 trials is plotted as solid lines, and the shading bounded by dotted lines denotes the standard deviation. By combining L=0 and L=2 components (in green), the standard deviation in axon diameter estimates decreases compared to using L=0 components alone (in orange), especially in the small diameter regime.

Figure 5. Estimated apparent axon diameter map of in vivo human data. A sagittal slice through the right corticospinal tract was shown for the map calculated using L=0 components only (left) and that obtained by combining L=0 and L=2 components (right). Overall, the two maps show similar patterns, but there is a subtle improvement in contrast between the larger apparent axon diameter in the corticospinal tract and surrounding white matter on the L0+L2 map, indicating that including higher-order spherical harmonics may provide better sensitivity to axon diameter compared to the spherical mean approach.

Proc. Intl. Soc. Mag. Reson. Med. 27 (2019)