Comparative analysis of signal models for microscopic fractional anisotropy estimation using q-space trajectory encoding

Microscopic diﬀusion anisotropy imaging using diﬀusion-weighted MRI and multidimensional diﬀusion encoding is a promising method for quantifying clinically and scientiﬁcally relevant microstructural properties of neural tissue. Several methods for estimating microscopic fractional anisotropy (μFA), a normalized measure of microscopic diﬀusion anisotropy, have been introduced but the diﬀerences between the methods have received little attention thus far. In this study, the accuracy and precision of μFA estimation using q-space trajectory encoding and diﬀerent signal models were assessed using imaging experiments and simulations. Three healthy volunteers and a microﬁbre phantom were imaged with ﬁve non-zero b-values and gradient waveforms encoding linear and spherical b-tensors. Since the ground-truth μFA was unknown in the imaging experiments, Monte Carlo random walk simulations were performed using axon-mimicking ﬁbres for which the ground truth was known. Furthermore, parameter bias due to time-dependent diﬀusion was quantiﬁed by repeating the simulations with tuned waveforms, which have similar power spectra, and with triple diﬀusion encoding, which, unlike q-space trajectory encoding, is not based on the assumption of time-independent diﬀusion. The truncated cumulant expansion of the powder-averaged signal, gamma-distributed diﬀusivities assumption, and q-space trajectory imaging, a generalization of the truncated cumulant expansion to individual signals, were used to estimate μFA. The gamma-distributed diﬀusivities assumption consistently resulted in greater μFA values than the second order cumulant expansion, 0.1 greater when averaged over the whole brain. In the simulations, the generalized cumulant expansion provided the most accurate estimates. Importantly, although time-dependent diﬀusion caused signiﬁcant overestimation of μFA using all the studied methods, the simulations suggest that the resulting bias in μFA is less than 0.1 in human white matter.


Introduction
Diffusion-weighted magnetic resonance imaging (dMRI) has become firmly established as the MRI technique of choice for quantifyhttps://doi.org/10.1016/j.neuroimage.2021.ing neural tissue's microstructural properties in vivo ( Johansen-Berg and Behrens, 2013 ).Since the pioneering work by Basser et al. (1994) , diffusion tensor imaging (DTI) parameters quantifying the magnitude and shape of the voxel-level diffusion tensor have been widely applied in characterizing brain development and pathologies ( Assaf et al., 2019 ).Despite its utility, DTI suffers from various well-recognized limitations that have been addressed by the subsequent advances in the field.Because orientation dispersion of anisotropic neurites can result in an isotropic diffusion tensor, several methods have been developed to obtain a more accurate estimate of the fibre orientation distribution, e.g., ( Tournier et al., 2007;Tuch, 2004;Wedeen et al., 2008 ).Another major limitation of DTI is that it does not provide a good fit to data in experiments with moderate to high diffusion-weighting (roughly  < ms ∕μm 2 in the brain) when the signal as a function of b-value clearly deviates from a monoexponential decay, revealing that the voxel-level diffusion propagator is not Gaussian, especially in white matter ( Jensen et al., 2005 ).More complicated signal models capturing the deviation from the monoexponential decay, e.g., ( Clark et al., 2002;Jensen et al., 2005;Novikov and Kiselev, 2010 ), as well as a plenitude of microstructural models, e.g., ( Jespersen et al., 2007;Sotiropoulos et al., 2012;Zhang et al., 2012 ), have been developed to obtain a better fit to data and additional information on tissue microstructure.Comprehensive overviews of the state of the field are provided by Jelescu and Budde (2017) and Novikov et al. (2019) , for example.
The methods referred to above belong to a class of single diffusion encoding (SDE) methods ( Shemesh et al., 2016 ) as they are based on acquisitions measuring the displacements of water molecules along a single dimension, i.e., the direction of the diffusion encoding magnetic field gradient.Despite being sensitive to several microstructural changes, SDE methods are fundamentally limited because they confound two major sources of voxel-level non-Gaussian diffusion, namely, the orientation dispersion of anisotropic diffusion and the size variance of microscopic diffusion environments, resulting in a lack of specificity ( Westin et al., 2016 ).This fundamental degeneracy of SDE acquisitions also prevents accurate measurement of microscopic diffusion anisotropy without substantial a priori information about tissue microstructure ( Henriques et al., 2019;Ianu ş et al., 2016;Lampinen et al., 2017 ).Multidimensional diffusion encoding (MDE), on the other hand, renders the dMRI signal sensitive to the displacements of water molecules along two or three dimensions.Over recent years, MDE methods, such as double diffusion encoding (DDE) and q-space trajectory encoding (QTE), have gained significant attention for their ability to resolve the fundamental degeneracy in data acquired with conventional SDE methods, capturing clinically and scientifically relevant information about tissue microstructure ( Cory et al., 1990;Eriksson et al., 2013;Henriques et al., 2020b;Jespersen et al., 2013;Lasi č et al., 2014;Mitra, 1995;Topgaard, 2017;Westin et al., 2016 ).
DDE waveforms consist of two trapezoidal pulse pairs separated by a mixing time and therefore require rather long echo times, resulting in a diminished signal.Using clinical whole-body scanners, asymmetric DDE waveforms can also produce significant concomitant gradients in spin echo experiments, causing artefacts and signal dropout ( Szczepankiewicz et al., 2019c ), an issue that can be avoided by using double spin echo ( Koch and Finsterbusch, 2008 ).QTE waveforms are designed to enable spin echo experiments with shorter echo times by following an optimized trajectory in q-space to define a b-tensor with the desired shape, orientation, and magnitude ( Sjölund et al., 2015 ).This makes it feasible to perform MDE experiments using more limited scanner hardware ( Szczepankiewicz et al., 2019b ), albeit at the cost of assuming diffusion time dependency to be negligible because QTE waveforms lack a well-defined scalar diffusion time.If this assumption is not valid, a discrepancy between the diffusion times of the gradient waveforms encoding different b-tensors shapes results in biased parameter estimates, as shown by Lundell et al. (2019Lundell et al. ( , 2017) ) , who proposed reducing the bias by using so-called tuned waveforms with similar power spectra at least along one direction.Furthermore, isotropic dif-fusion encoding with QTE may be orientationally variant due to significant time-dependent microscopic diffusion anisotropy ( Jespersen et al., 2019;de Swiet and Mitra, 1996 ).
By combining MDE acquisitions with different b-tensor shapes, it is possible to estimate microscopic diffusion anisotropy, i.e., the average anisotropy of the microscopic diffusion environments, referred to as compartments hereafter, irrespective of their orientation dispersion.Since microscopic diffusion anisotropy does not depend on the orientation dispersion of axons, it may be a clinically valuable biomarker for axonal degeneration.Indeed, promising results of the application of microscopic fractional anisotropy (μFA) ( Lasi č et al., 2014 ), a normalized measure of microscopic diffusion anisotropy that is equivalent to conventional fractional anisotropy (FA) if the compartments are aligned, have been reported in imaging brain tumours ( Szczepankiewicz et al., 2015;2016 ), multiple sclerosis lesions ( Andersen et al., 2020;Yang et al., 2018 ), and microstructural properties of white matter associated with schizophrenia ( Westin et al., 2016 ), epilepsy ( Lampinen et al., 2020 ), and aging and Parkinson's disease ( Kamiya et al., 2020 ).However, different studies have applied different methods for estimating μFA and the differences between the methods have received little attention thus far.Understanding how different μFA estimation methods relate to each other is crucial in planning future studies and interpreting previously reported results.Furthermore, reaching a consensus on the optimal method for measuring μFA is an essential requirement for future clinical studies.
In this study, imaging experiments and simulations were performed to compare the accuracy and precision of μFA estimation using the following signal models: truncated cumulant expansion of the powderaveraged signal ( Ianu ş et al., 2018;Jespersen et al., 2013;Lasi č et al., 2014;Nery et al., 2019;Yang et al., 2018 ), gamma-distributed apparent diffusivities ( Andersen et al., 2020;Lasi č et al., 2014;Szczepankiewicz et al., 2015;2016 ), and q-space trajectory imaging (QTI), a generalization of the truncated cumulant expansion to individual acquisitions using diffusion tensor distribution (DTD) as a microstructural model, ( Kamiya et al., 2020;Lampinen et al., 2020;Martin et al., 2020;Westin et al., 2016 ).The μFA maps in healthy volunteers calculated from the same data using the different methods were compared.The precision of the estimates was quantified by repeatedly imaging a microfibre phantom.Furthermore, since the ground-truth value of μFA was unknown in the imaging experiments, Monte Carlo random walk simulations were performed using axon-mimicking fibres for which the ground truth was known.The bias in μFA due to time-dependent diffusion was quantified by repeating the simulations with tuned waveforms and with triple diffusion encoding (TDE).

Theory
Conventional DTI represents diffusion at the voxel level by a symmetric positive definite 3 × 3 tensor  that captures Gaussian diffusion ( Basser et al., 1994 ).As the magnitude of diffusion-weighting is increased, the signal as a function of b-value starts to deviate from a monoexponential decay, revealing that the voxel-level diffusion propagator is not Gaussian, especially in white matter ( Jensen et al., 2005 ).At the voxel level, non-Gaussian diffusion is caused by intra-compartmental restricted diffusion and variance in diffusion properties across the compartments ( Henriques et al., 2020a ).The theory presented in this section is based on the assumption that the effects of intra-compartmental restricted diffusion are negligible.
If inter-compartmental exchange of water and intra-compartmental diffusion time dependence are negligible, the tissue in a voxel can be represented by a distribution of microscopic diffusion tensors representing the compartments ( Westin et al., 2014 ).In this case, the diffusionweighted signal can be expressed as where  0 is the signal without diffusion-weighting,   is a microscopic diffusion tensor,  6D is the normalized distribution of microscopic diffusion tensors in the voxel,  is the b-tensor used in the acquisition, ∶ denotes the generalized scalar product (  ∶  = ∑ 3  =1 ∑ 3 =1     ), and integration is performed over all symmetric positive definite 3 × 3 tensors ( Westin et al., 2016 ).QTE is based on the idea that if intracompartmental diffusion is Gaussian, the diffusion-weighted signal attenuation depends only on  , assuming that all the other imaging parameters are kept constant ( Westin et al., 2016 ).
Since recovering the six-dimensional microscopic diffusion tensor distribution is an ill-posed problem ( Topgaard, 2017 ), several methods have been developed for estimating relevant properties of it.One of these is μFA, a normalized measure of the average eigenvalue variance of the microscopic diffusion tensors ( Jespersen et al., 2013;Lasi č et al., 2014;Szczepankiewicz et al., 2016 ): where Var  ( ) denotes an operator that calculates the eigenvalue variance of a tensor, Tr ( ) denotes the trace of a tensor, and ⟨ ⟩ denotes averaging over the DTD.

Powder-averaged signal
The powder-averaged signal, i.e., the signal averaged over all possible diffusion encoding directions, is orientationally invariant and can be expressed as where  is the distribution of apparent microscopic diffusivities in the voxel,  is the b-value used in the acquisition, and   is the apparent microscopic diffusivity ( Lasi č et al., 2014 ). depends on the shape of the b-tensor used in the acquisition.The powder-averaged signal can be accurately estimated by averaging over the acquired diffusion encoding directions, given that the signal-to-noise ratio (SNR) is high enough and a sufficient number of directions is used ( Jespersen et al., 2013;Szczepankiewicz et al., 2019b ).

Second order cumulant expansion
Using the cumulant expansion up to the second order in  , the powder-averaged signal can be expressed as where MD = ⟨Tr (   )∕3 ⟩ is the voxel-level mean diffusivity and  is the second central moment, i.e., variance, of  ( Lasi č et al., 2014 ).By combining acquisitions with different b-tensor shapes, the total variance in apparent microscopic diffusivities can be decomposed into anisotropic and isotropic variances,  aniso and  iso , respectively ( Szczepankiewicz et al., 2016 ). aniso is proportional to the average eigenvalue variance of the microscopic diffusion tensors and  iso is equal to the variance of mean diffusivities of the microscopic diffusion tensors: In experiments with linear and spherical b-tensors, such as in this study, the total variance relates to anisotropic and isotropic variances as where the subscripts LTE and STE represent acquisitions with linear tensor encoding and spherical tensor encoding, respectively ( Szczepankiewicz et al., 2016 ).Therefore, μFA can be estimated by fitting a signal model, e.g, Eq. ( 4) , to the estimated powder-averaged signals acquired with linear and spherical b-tensors.

Higher-order correction
The limitation of estimating μFA as described above is that Eq. ( 4) is strictly valid as  → 0 .However, the applied b-value has to be sufficiently high for voxel-level non-Gaussian diffusion to be measurable.The accuracy of μFA estimation using the truncated cumulant expansion and DDE has been shown to improve if an additional term is introduced to account for higher-order effects ( Ianu ş et al., 2018 ): where  3 ∈ ℝ .This correction, which is applied to the signal model representing acquisitions with linear encoding, was introduced to correct for higher-order effects in the difference between the signals acquired with different b-tensor shapes.It can reduce the b-value-dependent underestimation of μFA ( Ianu ş et al., 2018 ).

Gamma-distributed apparent diffusivities
The distribution of apparent microscopic diffusivities can also be a priori assumed to be such that enables Eq. ( 3) to be analytically evaluated, e.g., gamma-distributed ( Jensen and Helpern, 2010 ): where  ∈ ℝ + and  ∈ ℝ + are the shape and scale parameters of the gamma distribution, respectively, and Γ is the gamma function.Under this assumption, the powder-averaged signal can be expressed as where MD =  and  =  2 ( Lasi č et al., 2014 ).However, the resulting parameter estimates can be inaccurate if the distribution of apparent microscopic diffusivities is not gamma-distributed.

Q-space trajectory imaging
The second order cumulant expansion has also been generalized to describe individual acquisitions instead of the powder-averaged signal.In QTI ( Westin et al., 2016 ), the signal in a single acquisition is expressed as where ⊗ denotes the tensor outer product and ℂ is a 3 × 3 × 3 × 3 covariance tensor defined as The covariance tensor contains 21 unique elements and enables the estimation of several relevant properties of the DTD.μFA can be calculated as where  iso and  shear are the isotropic and shear tensors of rank 2 and shape 6 × 6 as described by Westin et al. (2016) .

μFA estimation
Four signal models for μFA estimation were compared: the cumulant expansion of the powder-averaged signal up to the second order in  ( Eq. ( 4) ) and with a higher-order correction ( Eq. ( 9) ), the gammadistributed diffusivities assumption ( Eq. ( 11) ), and QTI ( Eq. ( 12) ).For brevity, the methods will be referred to hereafter as the cumulant model, higher-order model, gamma model, and QTI.
The powder-averaged signal was estimated by averaging over the acquired diffusion encoding directions.The cumulant model, higher-order model, and gamma model were fit to the estimated powder-averaged signal using a non-linear least squares trust region reflective algorithm ( Branch et al., 1999 ) in Scipy ( Virtanen et al., 2020 ). 0 , MD,  aniso ,  iso , and  3 were used as fit parameters from which μFA was estimated using Eqs.( 2) , ( 5) , and (6) .The fit parameters were constrained to be non-negative real numbers.The initial values of the fit parameters were:  0 = average signal over acquisitions without diffusion-weighting, MD = 1 μm 2 ∕ ms ,  aniso = 0 . 1 μm 4 ∕ ms 2 ,  iso = 0 . 1 μm 4 ∕ ms 2 , and  3 = 0 .
QTI was fit to data using the following linear equation: where where  is the number of acquisitions and  , ℂ ,   , and (   ⊗   ) are column vectors in Voigt notation ( Westin et al., 2016 ).The diagonal matrix  with elements   =   was used to correct for heteroscedasticity in the log-transformed data ( Jones and Cercignani, 2010 ).The matrix inversion in Eq. ( 15) was performed using the Moore-Penrose pseudoinversion ( Strang and Borre, 1997 ) in Numpy ( Walt et al., 2011 ).μFA was calculated according to Eqs. ( 13) and ( 14) .

Imaging experiments
Data was acquired using a prototype spin echo sequence ( Szczepankiewicz et al., 2019a ) on a Siemens Magnetom Prisma 3T (Siemens Healthcare, Erlangen, Germany) with a maximum gradient strength of 80 mT/m, maximum slew rate of 200 T/m/s, and 64-channel head coil at Great Ormond Street Hospital, London, United Kingdom.Data was preprocessed with random matrix denoising ( Veraart et al., 2016 ) and Gibbs ringing correction ( Kellner et al., 2016 ) using MRtrix3 ( Tournier et al., 2019 ) and distortion correction using FSL's topup and eddy ( Andersson and Sotiropoulos, 2016;Jenkinson et al., 2012 ).
Twelve diffusion encoding directions were used for b-values 0.1, 0.5, and 1 ms/μm 2 and 32 directions for b-values 1.5 and 2 ms/μm 2 .The directions were distributed uniformly around the surface of a sphere by combining the vertex coordinates of the icosahedron and dodecahedron ( Westin et al., 2016 ).The protocol also included 15 images without diffusion-weighting, one of which had the phase encoding direction reversed.Other relevant imaging parameters were: voxel size = 2 × 2 × 2 mm 3 , FOV = 256 × 256 mm 2 , TE = 103 ms, and phase partial Fourier = 6/8.Data was acquired using numerically optimized ( Sjölund et al., 2015 ) and Maxwell-compensated ( Szczepankiewicz et al., 2019c ) gradient waveforms encoding linear and spherical b-tensors ( Fig. 1 A-B).To avoid peripheral nerve stimulation and gradient coil heating, the slew rate was constrained to a maximum of 65 T/m/s when calculating the gradient waveforms using the software provided by Sjölund et al. (2015) .The waveforms for both linear and spherical encoding were rotated according to the diffusion encoding directions ( Szczepankiewicz et al., 2019a ).
Under the Gaussian phase approximation, the diffusion time of a QTE gradient waveform can be quantified by calculating the power spectrum of the corresponding wave vector ( Stepi š nik, 1993 ): where PSD stands for power spectral density,  is the wave vector magnitude along a single direction, and  is the diffusion frequency.The power spectra of the Cartesian components of the unrotated waveforms used in the imaging experiments are shown in Fig. 1 E-F.Since the STE waveform exhibits spectral anisotropy, which may lead to an orientationally variant signal ( Jespersen et al., 2019 ), the orientational variance of the STE acquisition was quantified using the correlation between the FA values estimated from LTE and STE data.If the STE acquisition is orientationally invariant, non-zero values of FA calculated from the STE data are caused by noise and they are uncorrelated with the FA values calculated from the LTE data.DTI was fit separately to LTE and STE data using a weighted linear least squares fit in Dipy ( Garyfallidis et al., 2014 ).

Volunteer experiments
The experiments were approved by the UCL Research Ethics Committee.Each participants gave written and informed consent prior to the scan.
Three healthy adult volunteers (two females and one male with ages ranging from 27 to 42 years) were scanned with whole-brain coverage and the following imaging parameters: 60 slices, TR = 10 s, and 2 repetitions.The total scan time was approximately 40 minutes per volunteer.Areas of cerebrospinal fluid were excluded from the analysis by excluding the voxels where MD > 2.5 μm 2 /ms.
The agreement between the μFA maps was quantified using the concordance correlation coefficient   ( Lin, 1989 ), which quantifies the agreement between two variables measuring the same quantity.The values of   range from -1 to 1, with perfect agreement at 1.

Phantom experiments
A phantom consisting of highly hydrophilic hollow polycaprolactone microfibres ( Zhou et al., 2018 ) was imaged with the following imaging parameters: TR = 3 s, 10 slices, and 1 repetition.The phantom contained three regions of interest (ROI) with different fibre configurations: parallel fibres, perpendicularly crossing fibres, and fibres with random orientations.The inner diameter of the fibres was 9.9 ± 1.2 μm (mean ± standard deviation) in the ROI containing parallel and crossing fibres, and 7.8 ± 0.5 μm in the ROI containing randomly oriented fibres.Scanning electron microscope images illustrating the microstructure of the phantom are shown in Fig. 2 .
The acquisition protocol was repeated 11 times to study the precision of μFA.The coefficient of variation (CV) was calculated at each voxel as CV = ( ∕ ) ⋅ 100% , where  is the standard deviation of μFA estimated using a given method over the repeated acquisitions and  is the mean μFA averaged over all acquisitions and methods.A single value of  was used to not penalize the methods that produce lower values of μFA.

Simulation experiments
GPU-accelerated Monte Carlo random walk simulations were performed using Disimpy ( Kerkelä et al., 2020 ) to compare the signal models in a scenario where the ground truth was known and to assess μFA bias due to time-dependent diffusion.The functionality of the simulator has been confirmed by comparing the simulated signals to analytical solutions using simple geometries, such as cylinders and spheres.Validation examples are available on the documentation ( https://disimpy.readthedocs.io).

Simulated microstructure
The simulated microstructure, shown in Fig. 3 A, consisted of 1,713,598 triangles generated using ConFiG numerical phantom generator ( Callaghan et al., 2020 ) to represent a white matter region with Fig. 1.The diffusion encoding gradient waveforms used in the imaging experiments (A-B) and simulations (A-D), and the corresponding power spectra (E-H).It is important to note that the power spectra of the unrotated triple diffusion encoding waveforms are identical along  ,  and  but not along arbitrary directions.

Fig. 2.
Representative scanning electron microscope images illustrating the microstructure of the parallel microfibres phantom (A) and the randomly oriented microfibres phantom (B).The phantom with crossing microfibres was made from the same sample as the one with parallel microfibres.

Fig. 3. (A)
The simulated voxel containing 381 axonmimicking fibres.(B) A histogram of the fibre radii (mean along the fibre length).The dashed line denotes the probability distribution from which the target radii where sampled.The fibre radii skew towards smaller radii than the target distribution due to the fact that as the fibres grow in ConFiG, they often have to shrink to fit into the available space.dispersed axons.Palombo et al. (2019) have shown that the triangulation of a simulated cell membrane has a small effect on the simulated signal.The target orientation distribution of the fibres was drawn from a Watson distribution with  = 2 ( Mardia and Jupp, 2009 ), target radius distribution was drawn from a gamma distribution with mean of 1 μm and standard deviation of 0.1 μm, and target fibre density was 70%.The resulting voxel contained 381 irregular and impermeable axonmimicking fibres at a fibre density of 60%.The fibre radii skew towards smaller radii than the target distribution, as seen in Fig. 3 B, due to the fact that as the fibres grow in ConFiG, they often have to shrink to fit into the available space.Periodic boundary conditions were used, i.e., the random walkers that left the voxel encountered repeating identical copies of the shown microstructure.To ensure that the ends of the fibres matched, cut fibres were extended with mirrored copies of themselves in the directions aligned with the main fibre direction, and the final voxel included 50% of the mirrored copies.In the plane perpendicular to the direction of the fibres, the fibres were grown in a way that ensures periodicity.The volume of the simulated voxel was 39 × 39 × 32 μm 3 .

Simulation details
Three simulation experiments were performed.Simulation 1 was performed using QTE and simulations 2 and 3 using TDE.The TDE waveforms encoding spherical b-tensors consist of three gradient pulse pairs with orthogonal directions and with identical diffusion encoding times ( ), diffusion times ( Δ), and gradient magnitudes.Similarly to DDE, the difference between the onset of the subsequent gradient pulses with orthogonal directions is referred to as mixing time ( ).
In each simulation, 3 ⋅ 10 6 random walks were generated using diffusivity of 2 μm 2 /ms. 5 ⋅ 10 3 random walkers were randomly placed inside each fibre and 1.095 ⋅ 10 6 in the extra-axonal space.Step length  was 0.31 μm and duration  was 8.1 μs.When a random walker collided with a surface, it was elastically reflected off the collision point in such a way that the walkers total path length during  was equal to .

Simulation 1
Simulation 1 was performed using the same waveforms as the imaging experiments ( Fig. 1 A-B).10 4 time steps were used.

Simulation 2
Following the approach proposed by Lundell et al. (2019Lundell et al. ( , 2017) ) , Simulation 2 was performed with the waveforms shown in Figs. 1 A and 1 C (  = 35 ms, Δ =  = 44 ms) to reduce the parameter bias due to spectral anisotropy of the spherical waveform and the discrepancy between the diffusion times of the waveforms encoding linear and spherical btensors.The power spectra of the Cartesian components of the unrotated TDE waveform ( Fig. 1 G) matches that of the linear waveform ( Fig. 1 E).It is important to note that the power spectrum of the unrotated TDE waveform is identical along  ,  , and  but not along arbitrary directions.31,720 time steps were used to keep the step length consistent with Simulation 1.

Simulation 3
To further reduce the effects of time-dependent diffusion, Simulation 3 was performed with TDE waveforms with short pulses, and long diffusion and mixing times (  = 3 ms, Δ =  = 100 ms).The long mixing time allows the random walkers to more fully sample their environments, minimizing the correlations between their positions during the different gradient pulse pairs.This is an important requirement for μFA estimation with DDE and TDE, which does not require the assumption of time-independent Gaussian intra-compartmental diffusion ( Jespersen et al., 2013;2019 ).The waveform encoding spherical b-tensors is shown in Fig. 1 D. The waveform encoding linear b-tensors was a TDE waveform with the same direction for the three pulse pairs and with the same pulse length, diffusion time, and mixing time as the spherical waveform.62,195 time steps were used to keep the step length consistent with Simulation 1.

Noise addition
To study the precision of μFA, SNR of the simulated signals was lowered to 25. Signals with Rician noise ( Gudbjartsson and Patz, 1995 ) were generated as where  is the simulated signal without noise and  and  are randomly sampled from a normal distribution with zero mean and standard deviation .Here, SNR refers to the signal without diffusion-weighting divided by .Noise addition was performed 10 4 times.
To confirm that the noise due to finite number of spins was much less than the added noise, Simulation 1 was repeated 5 times with different pseudorandom number generator seeds.SNR, estimated as the signal without diffusion-weighting divided by the mean standard deviation of the signals generated with different seeds, was 5,027.

Ground-truth μFA
The ground-truth μFA was calculated according to Eq. ( 2) by separately fitting diffusion tensors to the LTE signals from each compartment, i.e., the individual fibres and extra-axonal space.The diffusion tensors were estimated by fitting diffusion and kurtosis tensors ( Jensen et al., 2005 ) to data using a weighted linear least squares algorithm in Dipy ( Garyfallidis et al., 2014 ).To vary the ground-truth μFA, three signal fractions for the spins in the intra-axonal space (  intra ) were

Volunteer experiments
Representative axial slices of the μFA maps in one of the volunteers are shown in Fig. 4 .A visual inspection reveals that the higher-order model and gamma model result in greater μFA values than the cumulant model and QTI.The best agreement, as measured by the concordance correlation coefficient, was found between the cumulant model and QTI (   = 0 .91 ), followed by the cumulant model and gamma model (   = 0 .85 ), higher-order model and gamma model (   = 0 .85 ), gamma model and QTI (   = 0 .81 ), higher-order model and QTI (   = 0 .64 ), and cumulant model and higher-order model (   = 0 .64 ).
Fig. 5 shows representative axial slices of the difference maps in one of the volunteers.The cumulant model and QTI produced very similar maps.μFA maps estimated using the higher-order model exhibited smaller contrast between gray matter and white matter than the other methods, as can be seen particularly clearly in Fig. 5 D. The gamma model produced greater values of μFA than the cumulant model and QTI consistently across the brain.In white matter, the gamma model also produced greater values of μFA than the higher-order model.
The Bland-Altman plots in Fig. 6 show the voxel-wise differences between the methods' μFA estimates against their mean values across all volunteers.The good agreement between the cumulant model and QTI can be seen in the small mean difference (0.02) and the relatively small width of the 95% central range of the differences (0.18).Because  3 was constrained to be negative, the higher-order model resulted in higher values of μFA than the cumulant model in each voxel.On average, the higher-order model produced greater values of μFA than the other methods, especially at lower values of μFA, which are found in gray matter.When averaged over the whole brain, μFA values produced by the gamma model were 0.12 greater than the cumulant model and 0.10 greater than QTI.The whiskers extend to the furthest data point within 1.5 ⋅ IQR from the 1st or the 3rd quartile.
A significant correlation (  2 = 0 .13 ,  < 10 −3 ) was found between the FA maps (provided in the supplementary information) calculated from the data acquired with linear and spherical b-tensors, revealing that the acquisition with isotropic diffusion-weighting was moderately orientationally variant.
The average SNR, as quantified by the mean signal divided by the standard deviation over the images without diffusion-weighting, was 28.

Phantom experiments
μFA values averaged over the phantom ROIs from the repeated acquisitions are shown in Fig. 7 A-C.The gamma model consistently resulted in the greatest values of μFA, followed by the higher-order model, cumulant model, and QTI.The differences in mean μFA calculated using the different signal models were statistically significant in all pair-wise comparisons using an independent samples t -test (  < 10 −3 ).Despite the smaller fibre radius, μFA was the lowest in the ROI with randomly oriented fibres due to the lower packing density of the microfibres.
The distributions of the voxel-level CVs are shown in Fig. 7 D-F.The median CVs were well below 2%, indicating good repeatability in voxels with high μFA.The higher-order model stood out for having the lowest repeatability, especially in the phantom ROIs with crossing or randomly oriented fibres.A significant correlation (  2 = 0 .11 ,  < 10 −3 ) was found between the FA maps calculated from the data acquired with waveforms encoding linear and spherical b-tensors.
The average SNR was 123 as quantified by the mean signal divided by the standard deviation over the images without diffusion-weighting.

Simulation experiments
The results of the simulation experiments are presented in Fig. 8 .The ground-truth μFA was 0.34 for  intra = 0.2, 0.59 for  intra = 0.6, and 0.97 for  intra = 1.Furthermore, Bland-Altman plots of the simulation results are presented in Fig. 9 .
In Simulation 1, all methods overestimated μFA and the bias was greater at lower values of ground-truth μFA ( Fig. 8 A-C).In Simulation 2, the application of the TDE waveform whose power spectrum along three orthogonal directions matches that of the linear waveform reduced the overestimation of μFA for all methods ( Fig. 8 D-F).In Simulation 3, the application of the TDE waveforms with short pulses, long diffusion times, and long mixing times further reduced the bias ( Fig. 8 G-I).These results show that both the discrepancy between the diffusion times of the waveforms encoding different b-tensor shapes, and the correlations between the spins' positions during the diffusion encoding along different directions can result in overestimation of μFA.However, the simulations suggest that the bias is small in white matter, where μFA is high.The average positive bias caused by time-dependent diffusion, i.e., the difference between the expected μFA in simulations 1 and 3, was 0.11 for  intra = 0.2, 0.09 for  intra = 0.6, and 0.04 for  intra = 1.
The accuracy of the signal models for μFA estimation was assessed using the results from Simulation 3 which are not confounded by the bias induced by time-dependent diffusion.The mean squared error (MSE) between the expected value of μFA and the ground-truth value was 1.8 ⋅ 10 −3 for the cumulant model, 5.4 ⋅ 10 −3 for the higher-order model, 12.7 ⋅ 10 −3 for the gamma model and 1.6 ⋅ 10 −3 for QTI.In terms of precision, the average CV of μFA was 6.6% for the cumulant model, 10.6% for the higher-order model, 6.6% for the gamma model, and 8.5% for QTI.However, μFA is heteroscedastic and the precision of each method strongly depended on the ground-truth μFA.
The higher-order model is motivated by the b-value-dependent underestimation of μFA using the cumulant model.Indeed, the higherorder model produced more accurate estimates of μFA when the cumulant model underestimated μFA.However, because  3 was constrained to be non-negative, the higher-order model could not produce lower values of μFA and thus it provided even more inaccurate estimates in situations where the cumulant model overestimated μFA.Allowing  3 to take negative values improved the accuracy (MSE = 2.0 ⋅ 10 −3 .)but greatly destabilized the fit (mean CV = 24%).In most cases, at the bvalues used in the simulations, the underestimation due to the invalidity of the cumulant expansion of the powder-averaged signal was not strong enough to cancel the overestimation.

Discussion
Over recent years, MDE methods have gained significant attention for their ability to disentangle the orientation dispersion of anisotropic microscopic diffusion environments and the size variance of microscopic diffusion environments, enabling microscopic diffusion anisotropy to be measured without a priori information about tissue microstructure.Several microscopic diffusion anisotropy estimation methods have been introduced ( de Almeida Martins and Topgaard, 2016;Henriques et al., 2020a;Ianu ş et al., 2018;Jespersen et al., 2013;Lasi č et al., 2014;Westin et al., 2016 ), yet the differences between the methods have received little attention thus far.The purpose of this study was to assess the accuracy and precision of μFA estimates calculated using different signal models.μFA was chosen as the metric of interest for its connection to conventional FA that is well known in the neuroscience community.The cumulant model, gamma model, and QTI were chosen as the signal models of interest for their numerous applications in recent neuroimaging studies ( Andersen et al., 2020;Kamiya et al., 2020;Lampinen et al., 2020;Szczepankiewicz et al., 2015;2016;Westin et al., 2016;Yang et al., 2018 ).The higher-order model was included in the study for it being a simple extension to the cumulant model that has been shown to improve the accuracy of μFA estimation in animal studies using DDE ( Ianu ş et al., 2018 ).To our knowledge, only two studies have previously focused on comparing μFA estimation methods: Ianu ş et al. ( 2015) compared the cumulant model using DDE to the gamma model using TDE, and Reymbaut et al. ( 2020) compared the gamma model, QTI, and Monte Carlo inversion methods by explicitly modelling the tissue as a distribution of microscopic diffusion tensors using Eq. ( 1) .
In another relevant study, Martin et al. (2020) estimated the covariance tensor using QTI and investigated how μFA relates to other microscopic diffusion anisotropy indices in terms of contrast-to-noise ratio.In this study, four μFA estimation methods were applied on MRI data acquired using a clinical whole-body scanner and on synthetic data generated by simulating time-dependent diffusion.
In this study, the gamma model suffered from the greatest bias in the simulation experiments and consistently produced greater μFA values than the other methods in healthy white matter, fibre phantom, and simulations.The disagreement between the gamma model and the other methods was expected as the gamma model assumes the distribution of apparent microscopic diffusivities to follow a gamma distribution whereas the other methods are based on the truncated cumulant expansion.It is important to emphasize that although the assumption of gamma-distributed axon radii may be justified in real tissues ( Lee et al., 2019 ), the radii of our simulated axons did not follow a gamma distribution ( Fig. 3 B).Furthermore, gamma-distributed axon radii is not a sufficient condition for the apparent diffusivities to follow a gamma distribution.Reymbaut et al. (2020) have also reported that the gamma model produces biased estimates when the signal strongly deviates from a monoexponential decay, as is the case in white matter.However, they also reported that QTI results in biased estimates in the case of significant size variance of microsocpic diffusion environments.Therefore, given the microstructure-dependent bias, the results presented here can not be generalized to all situations.Nevertheless, our results showed a consistent disagreement between the gamma model and the second order cumulant expansion, which must be taken into account when interpreting reported results.
In terms of precision, the higher-order model stood out for being the most sensitive to noise.This is most likely explained by the third order polynomial fit's sensitivity to signal variation.Allowing  3 to take negative values further destabilized the fit.The higher-order model's low precision probably prevents it from being useful in most human neuroimaging studies.In the phantom experiments, the voxel-specific CV of μFA from all methods was low with the majority CV values being well below 2%, indicating very good repeatability in voxels with high microscopic diffusion anisotropy.However, since the phantom contained highly anisotropic microfibres, this result was expected based on previous reports of how the precision of μFA is high in voxels with high μFA but quickly diminishes with a decreasing value of μFA ( Ianu ş et al., 2016;Kerkelä et al., 2020;Lasi č et al., 2014;Szczepankiewicz et al., 2015;2019b ).
The μFA estimation methods used here assume intra-compartmental kurtosis to be negligible, which can result in biased estimates in the case of time-dependent diffusion ( Henriques et al., 2020a;2021 ).Our simulation experiments showed that the combination of trapezoidal linear gradient waveforms and isotropic encoding with QTE can lead to a significant overestimation of μFA, a result that has been discussed in detail by Lundell et al. (2019) .As a solution, they proposed designing the linear waveform so that its power spectrum matches that of the spherical waveform at least along one direction.Indeed, in the simulations, the accuracy of μFA estimates was improved when the power spectra were matched.However, the simulations also showed that this is not sufficient for accurate estimation of μFA.The accuracy of μFA further improved when a TDE waveform with short pulses, long diffusion time, and long mixing time was applied, eliminating the correlations between the positions of the random walkers during the applications of the gradient pulse pairs, an important assumption in the theory of μFA estimation with DDE ( Jespersen et al., 2013 ).However, the simulated TDE waveforms are completely unrealistic in today's human neuroimaging experiments because of their duration and gradient strength requirements.Considering that the bias in μFA induced by time-dependent diffusion was small when the ground-truth μFA was high, as is the case in white matter, and the effects of time-dependent diffusion are easier to observe in Monte Carlo random walk simulations than in the brain ( Gyori et al., 2020 ), these results support using QTE for μFA estimation in human white matter.
This study has several limitations.Most importantly, the accuracy and precision of μFA estimation depend on the details of the acquisition protocol which was kept constant here.For instance, the effect of b-values and diffusion encoding directions was outside the scope of this study.Also, the problem of designing optimal gradient waveforms for MDE experiments, discussed in detail by Szczepankiewicz et al. (2020) , was outside the scope of this study.Furthermore, our acquisitions with the lowest b-value may be affected by perfusion effects that may cause parameter bias.However, the focus of this study was to assess the performance of the signal models relative to each other when applied on the same data.It also must be mentioned that the microfibre phantom used for studying the precision of μFA contains fibres that are too large to be representative of white matter microstructure.In the simulations, the effects of orientation dispersion, distribution of axon radii, packing density, and membrane permeability were not assessed because the simulated geometry was kept constant across the simulations.Also, the pe-riodic boundary conditions introduced unrealistic periodicity into the simulated microstructure.Nevertheless, the simulated fibre geometry, generated using the algorithm by Callaghan et al. (2020) , resembles real axons and enabled the μFA estimates to be compared to a known ground truth.Finally, although Lee et al. (2021Lee et al. ( , 2020) ) have shown that a larger number of random walkers may be necessary to accurately calculate the moments of the displacement distribution, we only analyzed the simulated signal, which had converged, as shown by our results that are in good agreement with previously reported results ( Callaghan et al., 2020;Hall and Alexander, 2009;Rafael-Patino et al., 2020 ).Future studies should address the mentioned limitations to facilitate the optimal use of μFA.

Conclusion
The gamma-distributed diffusivities assumption produced greater values of μFA than the second order cumulant expansion.The simulations showed that matching the power spectrum of the linear waveform to the power spectrum of a component of the spherical waveform is not sufficient for accurate estimation of μFA.Nevertheless, the simulations suggest that the bias in μFA caused by time-dependent diffusion is small in human white matter.

Fig. 4 .
Fig. 4. Representative axial slices of the μFA maps in one of the volunteers calculated using the cumulant model (A), higher-order model (B), gamma model (C), and QTI (D).

Fig. 5 .
Fig. 5. Representative axial slices of the difference maps highlighting the differences between the methods' μFA estimates in one of the volunteers.

Fig. 6 .
Fig. 6.Voxel-wise comparison of the μFA maps.The differences are plotted against their mean value.The solid lines depict the mean difference.The dashed lines show the 2.5th and the 97.5th percentiles of the distribution of the differences.Colour represents data point density.

Fig. 8 .
Fig. 8.The distributions of the μFA estimates calculated from the simulated signals after 10 4 repetitions of noise addition to lower SNR to 25.   is the signal fraction of the intra-axonal space.(A-C) Results of Simulation 1. (D-F) Results of Simulation 2. (G-I) Results of Simulation 3. The dashed line denotes the ground-truth μFA.The blue line denotes the median.The box spans the interquartile range.The whiskers extend to the furthest data point within 1.5 ⋅ IQR from the 1st or the 3rd quartile.

Fig. 9 .
Fig. 9. Comparison of the μFA estimates calculated from the simulated data.The differences are plotted against their mean value.The solid lines depict the mean difference.The dashed lines show the 2.5th and the 97.5th percentiles of the distribution of the differences.Colour represents data point density.