Accuracy and precision of statistical descriptors obtained from multidimensional diffusion signal inversion algorithms

In biological tissues, typical MRI voxels comprise multiple microscopic environments, the local organization of which can be captured by microscopic diffusion tensors. The measured diffusion MRI signal can, therefore, be written as the multidimensional Laplace transform of an intravoxel diffusion tensor distribution (DTD). Tensor‐valued diffusion encoding schemes have been designed to probe specific features of the DTD, and several algorithms have been introduced to invert such data and estimate statistical descriptors of the DTD, such as the mean diffusivity, the variance of isotropic diffusivities, and the mean squared diffusion anisotropy. However, the accuracy and precision of these estimations have not been assessed systematically and compared across methods. In this article, we perform and compare such estimations in silico for a one‐dimensional Gamma fit, a generalized two‐term cumulant approach, and two‐dimensional and four‐dimensional Monte‐Carlo‐based inversion techniques, using a clinically feasible tensor‐valued acquisition scheme. In particular, we compare their performance at different signal‐to‐noise ratios (SNRs) for voxel contents varying in terms of the aforementioned statistical descriptors, orientational order, and fractions of isotropic and anisotropic components. We find that all inversion techniques share similar precision (except for a lower precision of the two‐dimensional Monte Carlo inversion) but differ in terms of accuracy. While the Gamma fit exhibits infinite‐SNR biases when the signal deviates strongly from monoexponentiality and is unaffected by orientational order, the generalized cumulant approach shows infinite‐SNR biases when this deviation originates from the variance in isotropic diffusivities or from the low orientational order of anisotropic diffusion components. The two‐dimensional Monte Carlo inversion shows remarkable accuracy in all systems studied, given that the acquisition scheme possesses enough directions to yield a rotationally invariant powder average. The four‐dimensional Monte Carlo inversion presents no infinite‐SNR bias, but suffers significantly from noise in the data, while preserving good contrast in most systems investigated.

translational motion of water molecules. Restrictions to the measured motion are then used to infer the geometrical configuration of the tissues inside a given voxel. Throughout the signal-encoding process, each water molecule has time to explore a micrometer-scale volume determined by the local diffusivity and the geometry of the cell membranes. The volume explored is only a minor fraction of the available space within the macroscopic imaging voxel, and the dMRI signal is thus a combination of signals arising from multiple microscopic environments with potentially different diffusion properties. [1][2][3][4] Classical diffusion observables such as mean diffusivity or fractional anisotropy 5 contain information about the voxel-averaged set of microscopic environments, and can be used to detect changes in the microscopic tissue structure. While voxel-averaged information is readily accessible, resolving the signal contributions from the various microscopic domains is still an open problem that prevents the unambiguous interpretation of in vivo dMRI data in terms of specific tissue properties. [6][7][8][9] In the context of dMRI, the measured signal  probes diffusion processes over a specific observational time-scale that depends on the choice of experimental time parameters. For a given observational time-scale, a common description of the subvoxel composition of heterogeneous tissues is obtained by considering the combined non-Gaussian diffusion effects of restriction and exchange and approximating the signal decay by a continuous weighted sum of exponential decays, 6,10-29 yielding the following multidimensional Laplace transform: where  0 is the signal acquired when the diffusion encoding gradients have zero amplitude, (D) is the distribution of apparent diffusion tensors D or ''diffusion tensor distribution'' (DTD) describing an apparent collection of subvoxel diffusion domains, 13 and b is the diffusion encoding tensor (or -tensor), an experimental variable that is completely defined by the diffusion-weighting magnetic field gradients. 22,30,31 The -tensor was introduced in earlier works as a -matrix accounting for cross terms between diffusion and imaging gradients. [32][33][34] Here, '':'' denotes the Frobenius inner product b ∶ D = ∑ i ∑ j ij ij , and integration spans the space of symmetric positive-semidefinite tensors. Changing the observational time-scale may very well lead to a different set of exponential decays as a result of restricted diffusion 35 and exchange, 36 which implies that the measured DTD may depend on the time-scales or spectral contents of the diffusion-encoding gradients. Even though such time-dependent effects have been measured in human brain white matter, [37][38][39][40][41][42] spinal cord, 43,44 and prostate, 45,46 using stimulated echo or oscillating gradient sequences specifically designed for varying the observational time-scale over extended ranges, the above DTD description holds for the limited range of time-scales probed by clinical dMRI experiments in the brain. 44,[47][48][49][50][51][52][53][54] While (D) provides a simple, and experimentally accessible, description of intravoxel tissue heterogeneity, the inversion of Equation (1) constitutes a challenging ill-posed problem, where clearly distinct distributions may result in signal decays that are indistinguishable within the experimental noise, thus complicating the detailed interpretation of dMRI data. The difficulties posed by the inversion of Equation (1) can be circumvented by fitting the acquired data to intricate signal models, which consider a predefined number of microscopic components and impose a priori constraints on the underlying diffusion properties. 19,[23][24][25]28,51,55,56 Alternatively, one can add complexity at the signal acquisition stage in order to encode more information into the diffusion signal and hopefully reduce the number of prior assumptions used to quantify the subvoxel heterogeneity. 31,42,50,57 This idea has motivated the development of tensor-valued diffusion-encoding techniques, 31,58 a number of multidimensional acquisition protocols wherein gradient waveforms are carefully designed to target specific features of (D).
The specificity of tensor-valued diffusion encoding is clarified by expanding the Frobenius inner product in Equation (1). To do so, let us assume a distribution of axisymmetric microscopic diffusion tensors D, the axial and radial diffusivities || and ⟂ of which, respectively, can be combined to define measures of isotropic diffusivity D iso = (D || + 2D ⟂ )∕3 and normalized diffusion anisotropy D Δ = (D || − D ⟂ )∕(3D iso ) ∈ [−0.5, 1]. 30 The iso and Δ parameters can be interpreted as measures of the ''size'' and ''shape'', respectively, of D. Also imposing axisymmetry on the diffusion where = Tr(b) is the usual -value, Δ ∈ [−0.5, 1] is the normalized anisotropy of b, and 2 ( ) = (3 2 − 1)∕2 is the second Legendre polynomial.
The angle defines the relative orientation between the diffusion tensor's main axis ( , ) and the encoding tensor's main axis (Θ, Φ), and is obtained from the spherical law of cosines cos = cos cos Θ + sin sin Θ cos( − Φ). Equation (2) shows clearly that the Δ parameter controls the effects of diffusion tensor anisotropy and orientation on the measured signal decay. For instance, setting Δ = 0 isolates the effects of iso , while data acquired at Δ = 1 contain the largest possible influence of diffusion anisotropy and orientation. Building upon the specificity of the Δ parameter, tensor-valued diffusion encoding techniques compare data acquired at arbitrary b ''shapes'' in order to resolve the effects of iso , Δ , and diffusion tensor orientation. By contrast, conventional diffusion encoding schemes are limited to a single Δ = 1 value and do not allow control over the influence of diffusion anisotropy on the acquired signal. Figure 1A presents in silico data from two different voxels that are isotropic at the macroscopic scale: one containing a heterogeneous isotropic system (top row), and another one comprising an anisotropic system with low orientational order (bottom row). The displayed data were simulated using a simple protocol comprising multiple -tensors b with varying ''size'' ( ), but constant ''shape'' ( Δ ) and orientation (Θ, Φ). Signals acquired with such protocols are naturally written as the Laplace transform of a 1D distribution of scalar effective diffusivities ( eff ): 30 Comparison in silico between conventional diffusion data and tensor-valued diffusion data in the context of Laplace inversions (SNR = 100). A, From left to right: voxel contents, acquired signal ∕ 0 for multiple -values and constant linear ''shape'' Δ = 1 and orientation (Θ, Φ) (black points), and theoretical distributions of diffusivities ( eff ) (black lines) associated with different isotropic components (top) and randomly oriented identical anisotropic components (bottom). The isotropic voxel content consists of two equiprobable populations, one with D iso = 0.2 μm 2 /ms and the other with iso = 1.6 μm 2 /ms. The anisotropic voxel content is made of components with D iso = 0.9 μm 2 /ms and Δ = 0.89. The results of a standard Laplace inversion algorithm 59,61 are shown with blue lines. B, Results of the same Laplace inversion algorithm 59,61 on orientationally averaged tensor-valued diffusion-encoded data ( , Δ ) Equation (7), inverted via Equation (8) and simulated for a voxel consisting of an equiprobable combination of the two cases in panel A. Such inversion yields a 2D DTD of isotropic diffusivities D iso and anisotropic ratios || ∕ ⟂ = (1 + 2 Δ )∕(1 − Δ ). 66 The signals are shown for Δ = 0, 0.5, and 1 (from bottom to top, see purple arrow). C, Comparison between the inversions of conventional and tensor-valued diffusion datasets for a voxel containing an isotropic component and an anisotropic component with low orientational order. The isotropic component (25% of the voxel content) has D iso = 2.7 μm 2 /ms and the anisotropic component (75% of the voxel content) has D iso = 0.9 μm 2 /ms and Δ = 0.9. The orientational order of the anisotropic component is set by a Watson distribution of order parameter OP = E[ 2 (cos )] = 0.2, where 2 ( ) = (3 2 − 1)∕2 is the second Legendre polynomial and is the angle between a component and the main eigenvector of the Saupe order tensor. 69,100 A fair comparison was ensured by considering the same number of input signals and the same SNR for both datasets. Inversions were performed using our 4D Monte Carlo inversion algorithm, 68 yielding a 4D DTD where the unit vectors giving the orientations of the microscopic components are color-coded via the [red, green, blue] ≡ [| sin cos |, | sin sin |, | cos |] convention. Ground truth is shown using black contours where eff ≡ iso [1 + 2 Δ Δ 2 (cos )]. Because the two voxels comprise clearly distinct microscopic structures, their theoretical ( eff ) distributions have markedly different functional forms. However, in spite of their distinct distributions ( eff ), the two intravoxel contents yield virtually indistinguishable Δ = 1 signals and thus cannot be differentiated with conventional diffusion encoding schemes. Moreover, Laplace inversion algorithms 59-61 cannot reproduce a complex broad distribution ( eff ) such as the one that characterizes the anisotropic system, and instead yield a solution that bears a higher resemblance to the theoretical distribution of the isotropic system. In the context of dMRI, the different peaks in ( eff ) may thus be wrongly interpreted as originating from intra-and extracellular water populations. Intuitively, these limitations can be understood via the fact that the one-dimensionality ( ) of this simple acquisition protocol does not match the two relevant dimensions ( iso , Δ ) of the physical problem at hand. Figure 1B shows that tensor-valued diffusion encoding circumvents this problem by adding the encoding ''shape'' Δ to the acquisition dimensionality. [62][63][64][65][66] The displayed data were simulated using the orientational average of the diffusion signal (see Equation (7)) and inverted via Equation (8), that is, the orientationally averaged equivalent of Equation (1). Finally, Figure 1C underlines the fact that conventional acquisition schemes sampling multiple directions may not capture the composition of a system that contains isotropic and anisotropic components, even at infinite signal-to-noise ratio (SNR). In other words, this figure demonstrates that the unique information encoded along the Δ dimension cannot be replaced by a simple increase in the sampling density of conventional experimental variables. Instead, substituting multiple Δ points with more -values or orientations (Θ, Φ) results in noisy distributions that do not enable a clear separation of the various subvoxel components.
Tensor-valued encoding protocols can be used to encode more information into the acquired signal, but do not correct for the intrinsic ill-posed character of Equation (1). Hence, estimating (D) from multidimensional dMRI data acquired at low SNR and comprising a limited number of ( , Δ , Θ, Φ) points remains a daunting task in itself. This has led to the development of multiple signal inversion approaches that are mostly distinguished by the amount of information they attempt to extract from (b). 31 These techniques include nonparametric approaches that attempt to estimate the underlying diffusion tensor distribution, [66][67][68] and parametric approaches that target statistical descriptors of (D) instead of the full distribution. 22,69 Simple and accessible descriptors are the mean isotropic diffusivity the normalized mean squared anisotropyẼ and the variance of isotropic diffusivities where aniso = iso Δ , and E[ · ] and V[ · ] are the voxel-scale mean and variance, respectively, built from the DTD. These descriptors can be related to quantities derived in previous works: E[ iso ] is identical to the mean diffusivity (MD);Ẽ[ 2 aniso ] carries similar information to the microscopic anisotropy index (MA), 70 the fractional eccentricity (FE), 71 the microscopic anisotropy (μA), 72,73 the normalized difference between second moments Δ̃2 and microscopic fractional anisotropy (μFA), 69,72 the anisotropic variance A and anisotropic mean kurtosis (MK A ), 74,75 and the microscopic anisotropy ; 22 V[ iso ] yields similar information to the isotropic second moment iso 2 , 69 the isotropic variance I and isotropic mean kurtosis (MK I ), 74,75 and the normalized isotropic variance MD . 22 Note that the conventional fractional anisotropy (FA) 5 is a convolution of anisotropy and orientational order, that is, it cannot tease apart the mean squared anisotropic diffusivity E[ 2 aniso ] and the order parameter OP = E[ 2 (cos )], as shown in previous work. 69 In fact, tensor-valued diffusion encoding is essential to capture subvoxel anisotropy, 76 reinforcing the conclusions drawn above from Figure 1C. Nevertheless, some caution is required when applying arbitrary gradient waveforms to systems exhibiting restricted diffusion, as directional discrepancies between the frequency contents of these waveforms may cause misestimations of diffusion metrics. 42,77,78 While such misestimations seem challenging to observe on a clinical scanner, 79 they can still be mitigated via ''tuning'' of the spectral contents. 42,57 Here, we compare techniques commonly used to invert tensor-valued dMRI data by assessing the accuracy and precision of their estimations . We focus on four techniques: • for 2D orientationally averaged signals ( , Δ ), the Gamma distribution fitting 69 (Gamma) and the 2D Monte Carlo inversion 66 (MC-2D); • for 4D signals ( , Δ , Θ, Φ), the covariance tensor approximation 22 (Cov) and the 4D Monte Carlo inversion 67,68 (MC-4D).
The comparison is performed using in silico datasets mimicking biologically plausible scenarios, probed with acquisition protocols and SNR levels that are compatible with a clinical setting. Although previous contributions have already reviewed the approaches listed above, 31,58 no rigorous assessment of their performances under clinical experimental conditions has been made so far. In section 2, we lay down the theory of each of the signal inversion techniques investigated. In section 3, we explain how the comparison between these techniques has been carried out. The results are presented and discussed in section 4. Finally, an outlook is given in section 5.

SIGNAL INVERSION APPROACHES
This section offers a pedagogical review of the signal inversion methods investigated in this work, namely the Gamma distribution fitting, 69 the covariance tensor approximation, 22 and the 2D 66 and 4D 67,68 Monte Carlo inversions. Readers with a good understanding of these methods may head directly to section 3.

Two-term cumulant expansion of the powder-averaged signal
The dimensionality of the problem stated in Equation (1) can be reduced by averaging data acquired at multiple b orientations: 80,81 Known as ''powder-averaging'', 30,82 this procedure results in a 2D signal that is insensitive to the orientation distribution of the subvoxel environments, and can be written as where our notation of the functional form of the effective distribution  Δ depends explicitly on Δ and implicitly on the subvoxel microstructure. General analytical expressions of the powder-averaged signal have been derived recently. 83 This signal corresponds to an integral transform of a joint distribution of isotropic and anisotropic diffusivities, ( iso , Δ ), that can be estimated from the data (see Figure 1B). 66 However, much like the original 4D problem, the direct inversion of Equation (8) is a challenging problem that has motivated a number of alternative approaches. For instance, one of the most straightforward ways of extracting information from ( , Δ ) consists of approximating Equation (8) with the mean = ∫ This expansion is only valid for → 0, which explains why it is commonly cut at the second order, as higher-order terms are negligible in that limit, yielding the two-term cumulant expansion By fitting the measured signal to Equation (10), using  0 , , and 2 as fitting variables, one can target statistical measures of  Δ ( ) directly without having to either retrieve the full distribution or define its explicit functional form. From the cumulant expansion, it is clear that signal deviations from monoexponentiality are described at low -values by the second central moment 2 . These deviations may come from intravoxel variance in isotropic diffusivities or from anisotropy, and their effects on the powder-averaged signal can be teased apart by the encoding anisotropy Δ , as 2 can be written as 30 Within this formalism, one has Even though the cumulant approach seems very direct for microstructural studies, it is inherently limited to small values of . Importantly, this means that adding more terms to Equation (10) (in the spirit of Equation (9)) will not significantly improve the adequacy of the cumulant expansion to describe systems where higher-order moments of  Δ ( ) are significant, such as ensembles of orientationally dispersed anisotropic components (see the bottom panel of Figure 1A).

Gamma distribution fitting (Gamma)
The problem of inverting ( , Δ ) (see Equation (7)) can instead be formulated as finding a physically plausible approximation of  Δ . The Gamma distribution function has been shown to provide a good approximation: 85 where is the unitless shape parameter, is a scale parameter with diffusivity units, and Γ denotes the Gamma function. This analytical form of  Δ is practical, because it is described by a minimal number of parameters and possesses a simple analytical Laplace transform: ( , , )∕ 0 = (1 + ) − . Using the fact that the Gamma distribution gives = and 2 = 2 , one obtains the analytical expression 69 the fit of which provides and 2 (b Δ ) for microstructural estimation of the metrics Equations (12), (13), and (14). It is important to mention that the Gamma distribution fitting implemented in this work is performed on the entire dataset and features a soft Heaviside signal-weighting function that limits the fit to the range where the signal attenuation does not reach a reasonable noise floor, as used in previous works. 74,86

Covariance tensor approximation (Cov)
The two-term cumulant expansion Equation (10) can be generalized to non powder-averaged signals as the covariance tensor approximation: 22 with the voxel-scale averaged diffusion tensor ⟨D⟩, the outer tensor product ⊗ such that b ⊗2 = b ⊗ b, and the 6×6 covariance tensor C = ⟨D ⊗2 ⟩ − ⟨D⟩ ⊗2 in Mandel notation. This approximation is equivalent to considering a normal distribution of diffusion tensors. However, the fact that this distribution allows for negative-definite tensors implies that negative values can be found for the positive quantities described in section 3.2. The aforementioned metrics Equations (12), (13), and (14) can then be expressed in terms of tensor inner and outer products involving the voxel-scale averaged diffusion tensor and the covariance tensor. Indeed, let us introduce the isotropic tensors E iso = I 3 ∕3 and E iso = I 6 ∕3, where I is the × identity matrix. Using the Mandel notation, in which a 3×3 symmetric tensor Λ can be written as the column vector ( xx yy zz √ 2 yz √ 2 xz √ 2 xy ) T , one builds the bulk and shear modulus tensors, by analogy with the stress tensor in mechanics, and obtains The list of fitting parameters includes  0 , the six independent elements of ⟨D⟩, and the 21 independent elements of C, which makes a total of 28 parameters to estimate within this approximation. Consequently, the covariance tensor approximation requires the acquisition of multiple signals over a wide range of -tensor sizes, shapes, and orientations in order to be reliable. The main advantage of the covariance tensor approximation is that it yields very fast inversions, as Equation (17) can be recast into a linear problem tractable via mere matrix pseudoinversion. However, this approach's limitation to small values of b ∶ ⟨D⟩ prevents it from accurately capturing the microstructural features of systems where higher-order moments of the underlying DTD are significant, such as ensembles of orientationally dispersed anisotropic components (see the bottom panel of Figure 1A). Figure 2 shows that the Gamma (Equation (16)) and covariance (Equation (17)) approaches may not be able to capture the features of the DTD (D) accurately, even at infinite SNR. Used as forward models, they can yield biased signals; used to solve an inverse problem, they can yield biased statistical descriptors.

Monte Carlo inversion (MC-2D and MC-4D)
To facilitate its numerical inversion, the continuous distribution (D) can be approximated by a discrete vector w containing the weights of distinct diffusion tensors. This allows us to rewrite Equation (1) as where  is the th signal amplitude measured with the encoding tensor b , and is the weight of the th point from the discretized DTD space. The analysis D and acquisition b spaces are linked by the generalized inversion kernel , the analytical form of which is determined by the experimental design. For 4D signals (b) , the inversion kernel is determined by Equation (2), whereas for 2D powder-averaged signals ( , Δ ) the kernel can be shown to be equal to 30 with erf( · ) denoting the error function. Since the number of experimental points, , is typically higher than the number of unknowns, , Equation (?) represents an overdetermined system. The estimation of w is thus formulated as a nonnegative linear least-squares (NNLS) problem: where || · || 2 denotes the Euclidian norm, w is the × 1 sought-for probability vector, s is the × 1 vector containing the signal amplitude measurements, and K is the × matrix approximation of the appropriate inversion kernel. FIGURE 2 Comparison at infinite SNR between typical results yielded by the Gamma distribution fitting (orange), the covariance tensor approximation (green), and the ground-truth signals/statistical descriptors (black) obtained for the same voxel content as in Figure 1C. In the signal panels, the colored solid lines denote the signal decays estimated from the various methods and the colored dashed lines correspond to their initial slopes, which are regulated by the estimated mean diffusivities. In the square DTD panels, colored points represent the logarithms of E[ iso ] and E[ || ∕ ⟂ ], and the horizontal line reflects the standard deviation of isotropic diffusivities. Since Δ is not directly accessible for The performances of Gamma and Cov are illustrated here for both the forward problem (top), where the ground-truth statistical descriptors are used to compute the signals predicted by the inversion methods, and the inverse problem (bottom), where the ground-truth signals are inverted to estimate the statistical descriptors. The signals are shown for Δ = 0, 0.5, and 1 (from bottom to top, see purple arrows), as in Figure 1B Inverting an ill-posed system of linear equations such as the one described by Equation (25) is a rather common task for a wide array of scientific disciplines. 87 Within the field of magnetic resonance, this class of problems is particularly relevant, as the analysis of relaxation or diffusion data can generally be formulated as an inversion of a discrete vector describing the amplitudes of multiple exponential decays. The ubiquity of the problem explains the abundance of methods that can be used to invert Equation (25), in both general 59,87 and magnetic resonance 61,[88][89][90][91] literature. In this work, we opted for iterative Monte Carlo algorithms that were specially designed to handle high-dimensional datasets: [66][67][68] MC-2D and MC-4D. While MC-2D is used to estimate ( iso , Δ ) from powder-averaged signals, MC-4D is used to estimate ( iso , Δ , , ) from non-powder-averaged signals.
The Monte Carlo inversion process was described in detail in previous works, 67,68 and is schematized in Figure 3. As a first step, the algorithm

Simulated acquisition schemes
We considered an abbreviated b acquisition scheme that can be readily implemented in a standard clinical MRI scanner as a 6-7 minute protocol comprising 30 axial slices. Similar protocols have already been used for in vivo studies of brain tissue microstructure. 95 The acquisition scheme

Comparison steps
To ensure a fair comparison of the different inversion techniques, all readily available as Matlab code on GitHub, 93,94 we follow the steps below. 1. A system of interest is simulated by generating a set of ground-truth features { || , ⟂ , , , }, from which the ground-truth descriptors of interest and the ground-truth set of signals { gt, } are computed using the general signal decay Equation (1). This computation is in agreement with the conventional procedure for testing 1D or 2D Laplace inversion algorithms, where the ground-truth signal is calculated with the same kernel as the inversion and it is well known that infinite SNR is not by itself sufficient for recovering the input distribution. 59,61,89,91 2. Each signal inversion technique is run on an identical set of signals with added Rician noise: where and ′ denote random numbers drawn from a normal distribution with zero mean.

Step 2 is repeated 100 times to build up statistics on parameter estimation of E[ iso ],Ẽ[ 2 aniso ], and V[ iso ].
In this work, we investigated a clinically relevant SNR of 30 and the ideal infinite SNR. Since the Rician bias has been shown to be relevant only when the SNR approaches unity, 96 affecting the estimation of diffusion metrics, 97-99 the results presented in section 4 are identical to those obtained from signals with added Gaussian noise. While accuracy is quantified by the difference between the median and the ground-truth (bias), precision is quantified by the interquartile range (IQR). A ''good'' estimation is defined as one that presents both high accuracy and high precision.

4.1
Varying mean diffusivity Figure 5 shows the accuracy and precision of the investigated signal inversion algorithms for isotropic unimodal systems with various mean diffusivities E[ iso ], probed with the acquisition scheme described in section 3.  Figure 5 and all remaining figures of this work, a general trend is that any bias in estimating the mean diffusivity E[ iso ] is reflected directly as a bias in estimating the other isotropic descriptor, that is, the isotropic variance V[ iso ]. In contrast, the estimation of the mean normalized anisotropyẼ[ 2 aniso ] is quite independent of the other descriptor estimations. Figure 6 shows the accuracy and precision of the signal inversion algorithms investigated for isotropic bimodal systems that change in terms of isotropic  (with Cov allowing negative anisotropy parameters, as in Figure 5). The Monte Carlo inversions do not show any infinite-SNR bias, but do exhibit finite-SNR biases at low E[ iso ] that increase with increasing V[ iso ]. This can be understood by the fact that generating isotropic variance implies giving more weight to intravoxel diffusional components with isotropic diffusivities both lower and higher than E[ iso ]. Hence, at low E[ iso ] and high V[ iso ], the in silico dataset comprises slow-diffusing components that present little signal contrast with the moderate -values used in this study. This means that, at low SNR, the MC-4D method cannot tease apart the slow-diffusing components and overestimates E[ iso ], which in turn leads to an overestimation of both  Figure 7 shows the accuracy and precision of the signal inversion algorithms investigated for anisotropic unimodal systems with varying normalized mean anisotropyẼ[ 2 aniso ]; while panel A focuses on orientationally ordered systems, panel B presents orientationally dispersed systems. Notice that, unlike other measures of anisotropy such as fractional anisotropy (FA), 5Ẽ [ 2 aniso ] is not affected by the degree of orientational order and is directly proportional to the underlying subvoxel anisotropy. 69 All diffusional components in Figure 7 share the same value of iso , which, according to Equation (5) For the most anisotropic system, the first-order signal deviation from monoexponentiality is quantified as (4∕5)Ẽ[D 2 aniso ] E[D iso ] 2 = 0.29 μm 4 /ms 2 (see Equation (11)). Gamma exhibits similar results for both oriented and orientationally dispersed systems: it can quantify anisotropy accurately, except for extremely high values ofẼ[ 2 aniso ], where infinite-SNR biases are observed. Similar infinite-SNR biases were also observed for isotropic systems FIGURE 6 Performance of different inversion techniques in the characterization of two isotropic toy systems consisting of bimodal distributions of D iso with a constant mean isotropic diffusivity E[ iso ] of (A) 0.8 μm 2 /ms and (B) 2 μm 2 /ms. Each toy system was devised as a sum of two Gaussian distributions with identical standard deviations ( = 0.05 μm 2 /ms). Symbol/color conventions are identical to those of Figure 5 with low E[ iso ] and high V[ iso ] (see Figure 6A), which indicates that the Gamma model cannot capture large deviations from monoexponentiality accurately. Gamma is not affected by orientational dispersion, because the powder-averaged signal associated with the acquisition scheme described in section 3.1 yields sufficient rotational invariance. Cov shows the same trends observed in the Gamma analysis, except for its infinite-SNR biases in orientationally dispersed systems containing microscopic domains of sufficiently high anisotropy. This can be understood by the fact that such systems are characterized by a ''ski-slope'' distribution (see the bottom panel of Figure 1A), wherein higher-than-second-order moments, not captured accurately by a cumulant expansion, are nonnegligible. Regarding MC-2D, it shows satisfying accuracy for all descriptors FIGURE 7 Descriptors of interest estimated by the inversion algorithms investigated as a function of the normalized mean anisotropyẼ[ 2 aniso ] for anisotropic toy systems consisting of (A) aligned and (B) orientationally dispersed anisotropic components, each of which is described by a unimodal distribution of Δ (Gaussian distribution) with constant iso = 0.8 μm 2 /ms and constant standard-deviation-to-the-mean ratio of 0.1. Orientationally dispersed systems are characterized by the order parameter OP = E[ 2 (cos )] = E[(3cos 2 − 1)∕2] = 0.01. Notice that even though the first and third simulated systems exhibit the same normalized mean anisotropy, the former contains planar anisotropic components (subscript P) and the latter contains linear anisotropic components (subscript L). Symbol/color conventions are identical to those of Figure 5 but has low precision when estimating E[ iso ] and V[ iso ]. Besides this, it is worth mentioning that MC-2D is very sensitive to the rotational invariance of the powder-averaged signal, as severe infinite-SNR biases appear when the powder-averaged signal is not truly rotationally invariant. A slight infinite-SNR bias is found in Figure 7A, because the acquisition scheme described in section 3.1 does not possess a perfectly

4.4
Orientational order and crossings  100 and meningiomas. 74 Given the fact that the different systems studied share common E[ iso ] = 0.8 μm 2 /ms,Ẽ[ 2 aniso ] = 0.46, and V[ iso ] = 0, the system at OP = 1 matches a combination of the two last systems of Figure 7. Consequently, only changes occurring with decreasing orientational order remain to be discussed, as other trends were already discussed in section 4.3. Moreover, the effects of decreasing orientational order were observed to be analogous to those of increasing crossing angles between bundles of anisotropic components, meaning that the observations made in this section can also be used to understand the performance of the various analysis methods for systems consisting of crossing fibers.
The accuracy and precision of Gamma are unaffected by orientational order, as the powder average of signals associated with the acquisition scheme described in section 3.1 yields sufficient rotational invariance. Any bias of Gamma is here inherited from the high subvoxel anisotropy of the different systems, as shown in Figure 7. As for Cov, its infinite-SNR biases identified at low OP and highẼ[ 2 aniso ] in Figure 7B are now observed to develop progressively with decreasing OP, as the higher-than-second-order moments of the underlying DTD gain more relevance. MC-2D still presents the lowest precision, but exhibits an increased accuracy with decreasing orientational order that correlates with its infinite-SNR biases. This can be explained by the fact that the acquisition scheme described in section 3.1 does not possess a perfectly rotationally invariant powder average. This limitation becomes less and less relevant as the system itself becomes powder-like. Because the Gamma and MC-2D methods are both based on powder-averaged signals, one would expect the performance of both methods to be independent of OP. The fact that MC-2D shows a decreased accuracy at high OP, while Gamma does not, shows that the former method is more sensitive to errors in the powder-averaging procedure than the latter. Finally, MC-4D does not possess any infinite-SNR bias, but sees its finite-SNR accuracy decrease with decreasing orientational order. This is due to the fact that fitting becomes more challenging when the subvoxel microstructure comprises multiple relevant orientations. FIGURE 9 Descriptors of interest estimated by the inversion algorithms investigated as a function of the signal fraction f iso ∈ [0, 1] of an isotropic cerebrospinal-fluid-like structure (Gaussian distribution of D iso with mean 3 μm 2 /ms and standard deviation 0.1 μm 2 /ms), starting from a voxel containing a fiber-like structure. This anisotropic structure is made of a Watson distribution of anisotropic components with an internal distribution of || (Gaussian distribution centered at 1.77 μm 2 /ms) and ⟂ (Gaussian distribution centered at 0.31 μm 2 /ms). The diffusivities used come from multiple white-matter fiber estimations of the Magic DIAMOND model. 23,24,104 The Watson distribution is set to an intermediate order parameter OP = 0.4 to mimic some fiber dispersion. Symbol/color conventions are identical to those of Figure 5

4.5
Varying fractions of isotropic and anisotropic components Figure 9 presents the accuracy and precision of the signal inversion techniques investigated for systems shifting progressively from a single fiber-like voxel content to a cerebrospinal-fluid-like voxel content; these systems mimic the progression of damage in a single white matter system with intermediate orientational order. This shift is quantified by the signal fraction f iso of the isotropic environment. The changes occurring between these composite systems are summed up here: where E[ iso ] is expressed in μm 2 /ms and V[ iso ] is expressed in μm 4 /ms 2 . Going from system 1 to system 5, E[ iso ] increases and the mean anisotropy decreases. As for isotropic variance, it increases between systems 1 and 3 and decreases between systems 3 and 5. Here, all results related to accuracy and precision can be understood from previous sections, so that only the contrasts captured by the inversion techniques remain to be discussed.
All investigated methods capture changes in E[ iso ] as a function of iso accurately, with comparable levels of accuracy and precision. Similarly, all methods capture the monotonic decrease ofẼ[ 2 aniso ] with increasing iso . While they share comparable precisions, the highest finite-SNR accuracy is obtained by MC-2D and the lowest finite-SNR accuracy is yielded by MC-4D. Finally, only MC-2D and MC-4D capture the changes in V[ iso ] as a function of iso . They present, respectively, low and intermediate precision.

CONCLUSION AND OUTLOOK
We investigated how the accuracy and precision of multidimensional signal inversion algorithms behave in multiple systems varying in terms of mean diffusivity, variance of isotropic diffusivities, mean squared diffusion anisotropy, orientational order, and fractions of isotropic and anisotropic components. The Gamma distribution fitting (Gamma) quantifies changes induced by E[ iso ] and is robust to orientational dispersion, but does not capture large signal deviations from monoexponentiality, originating from either isotropic variance or anisotropy, accurately. It is important to mention that all observations made in this article depend critically on the design of the acquisition scheme. Indeed, besides the aforementioned effect of rotational invariance of the powder-averaged signal, preliminary work by the present authors indicates that the performance of the various methods is affected by the choice of acquired coordinates ( , Δ , Θ, Φ), especially for Monte Carlo inversion techniques.
For instance, removing the = 0 signal from the acquisition scheme chosen in the present study, which is not conventional in the clinical setting, causes MC-4D to exhibit infinite-SNR biases at high V[ iso ] in the systems shown in Figure 6. This is due to the fact that our two lowest non-zero -values, 100 and 700 s/mm 2 , cannot capture alone the initial signal decay at low -value. Interestingly, acquisition schemes containing planar tensors ( Δ = −0.5) have been shown to increase precision in parameter estimation in both biophysical models and signal representations of white-matter tissues. 101,102 Although outside the scope of this present work, a comprehensive study of the link between experimental design and inversion performance should be pursued nonetheless. Alternative analysis procedures and improvements to the methods studied should also be investigated.
First, developing a Gamma equivalent of the covariance tensor approximation, that is, considering a matrix-variate Gamma distribution 23,24,103,104 in order to constrain the solution space to physical solutions, could prove useful. Second, previous studies have shown that denoising techniques 105 can improve the precision of parameters estimated from the Laplace inversion of MRI relaxometry data. 106 Assuming that similar denoising algorithms can be translated directly to tensor-valued diffusion data, they may improve the performance of the algorithms studied in this article.
The MC-4D method in particular could gain a huge boost in accuracy, as it does not possess any infinite-SNR bias and most of its limitations stem from its sensitivity to noise in the data. Finally, regularization techniques have traditionally been used to constrain the solution space of various diffusion and relaxation signal inversion techniques. 88,90,107,108 Whether or not such regularization constraints improve the inversion of tensor-valued data is a question that remains to be addressed.
Despite the protocol dependence of our present study, let us propose the following take-home message. The Gamma distribution fitting grants a versatile inversion that retrieves accurate statistical descriptors of the intravoxel diffusion tensor distribution as long as the powder-averaged signal does not deviate severely from monoexponentiality. The covariance tensor approximation shares many limitations of the Gamma distribution fitting, but offers fast inversions that render it practical for immediate, yet preliminary, clinical assessment of tissue microstructure.
Nevertheless, one should avoid studying anisotropic tissues with low orientational order, such as edematous white-matter tissues and certain meningiomas. Regarding these two inversion methods, preliminary work suggests that their infinite-SNR biases are intrinsic, which indicates that only limited improvement will be achieved by optimizing the acquisition scheme. In contrast, Monte Carlo inversion techniques are slow, but do not present infinite-SNR biases when applied to reasonable clinical protocols, which suggests that their true potential is yet to be reached, upon reduction of the effect of noise. Nonetheless, the 2D Monte Carlo inversion, the accuracy of which depends strongly on the rotational invariance of the powder-averaged signal and the precision of which is remarkably lower than that of other inversion techniques, should be avoided. As for the 4D Monte Carlo inversion, it may still preserve contrast despite its finite-SNR biases and is the only one that can access fiber-specific information, 92 making it promising on the one hand for clinical assessment of tissue microstructure and on the other hand for future studies of both neurodevelopment and neurodegeneration.