A new framework for MR diffusion tensor distribution

The ability to characterize heterogeneous and anisotropic water diffusion processes within macroscopic MRI voxels non-invasively and in vivo is a desideratum in biology, neuroscience, and medicine. While an MRI voxel may contain approximately a microliter of tissue, our goal is to examine intravoxel diffusion processes on the order of picoliters. Here we propose a new theoretical framework and efficient experimental design to describe and measure such intravoxel structural heterogeneity and anisotropy. We assume that a constrained normal tensor-variate distribution (CNTVD) describes the variability of positive definite diffusion tensors within a voxel which extends its applicability to a wide range of b-values while preserving the richness of diffusion tensor distribution (DTD) paradigm unlike existing models. We introduce a new Monte Carlo (MC) scheme to synthesize realistic 6D DTD numerical phantoms and invert the MR signal. We show that the signal inversion is well-posed and estimate the CNTVD parameters parsimoniously by exploiting the different symmetries of the mean and covariance tensors of CNTVD. The robustness of the estimation pipeline is assessed by adding noise to calculated MR signals and compared with the ground truth. A family of invariant parameters and glyphs which characterize microscopic shape, size and orientation heterogeneity within a voxel are also presented.

1 Maximum entropy property and identifiability of the CNTVD Theorem 1. For a Borel-measurable set, C ⊆ R d , the constrained Gaussian distribution with density, where Z(µ, Ω) is the normalizing constant is the unique maximizer of the entropy H(p) = − C p(X) log p(X)dX among all probability densities supported by C with given first and second moments X k = E µ,Ω (X k |X ∈ C), X k X = E µ,Ω (X k X |X ∈ C), 1 ≤ k, ≤ d.
When C contains an open ball the parameters µ, Ω are uniquely determined by the first and second moments.
Proof. Consider the following Lagrangian with the constraint p(X) ≥ 0, where λ ∈ R, ν ∈ R d and Ξ ∈ R d×d symm are the Lagrange multipliers. A solution of the maximization problem corresponds necessarily to a zero of the variational derivatives, δL δp(X) = − log p(X) − 1 + λ + X · ν + XΞX = 0 ∀X ∈ C (3) which means that when the maximum entropy solution exists necessarily it has the form p(X) = exp λ − 1 + X · ν + XΞX 1(X ∈ C) (4) with λ corresponding to a normalizing constant and ν, Ξ such that the first and second moment constraints are satisfied. When a solution with the given moments exists, after reparametrization, the density is rewritten as When the support C is unbounded, necessarily Ω is positive definite on the recession cone of C, and when Ω is positive definite on R d , we obtain a constrained Gaussian distribution. When a maximum exists, the uniqueness of the maximizing parameters is shown below.
Suppose that p µ, Ω (X) is another constrained Gaussian distribution with support C satisfying the moment constraints and maximizing the entropy. Then since the distributions satisfy the same moment constraints, which implies which means that Since pμ ,Ω (X) > 0 on C, for some r ∈ R, which is a variety of dimension, (d − 1), when µ = µ or Ω = Ω, contradicting the assumption of C having full dimension 2 Relations between the diffusion tensor cumulants and the dis-

placement cumulants
Let X ∈ R d be a conditionally Gaussian displacement vector, with zero conditional mean and conditional covariance 2D, while D is a random symmetric positive definite matrix (diffusion tensor) with probability density p(D). The Characteristic Function (CF) of the Gaussian conditional law of X given D is given by while the unconditional CF is obtained by integrating w.r.t. p(D) the conditional CF and it coincides with the Laplace Transform (LT) of the diffusion tensor distribution p(D) at Q = q q.
The cumulant expansion of the CF of X is given by where all the odd moments and cumulants vanish since the law of X is reflection symmetric, while the cumulant expansion of the LT of D is given by , andβ is the restriction of the multi-index β to the off-diagonal entries. By comparing these expansions for rank-1 Q = q q, we see that where the sum in Equation (23) is over the pairings π of {j 1 , . . . , j 2n } into pairs π 1 , . . . , π n , in accordance with the law of total cumulance. The even cumulants of the displacement κ α (X) with |α| = 2n are, up to rescaling, symmetrizations of corresponding cumulants of the diffusion tensor κ β (D) with |β| = n .
By truncating the cumulant expansion Equation (20) of the LT of p(D) up to second order terms, we obtain the LT of the NTVD with mean D and Covariance(D ij , D k ) = C ijk , which can be used to model the signal decay in a multiple pulsed gradient experiment: For single pulsed gradient experiments Q = q q has rank 1, and Equation (24) coincides with the signal decay in Diffusion Kurtosis Imaging (DKI) [1], obtained by truncating the cumulant expansion of the CF of X (Equation (17)) up to 4th order terms, where are respectively the diffusion tensor and the diffusion kurtosis tensor. Thus the NTVD signal model extends the DKI signal model from single to multiple pulsed experiments, and up to rescaling, the diffusion kurtosis tensor is the symmetrization of the partially symmetric NTVD covariance tensor.
Both NTVD and DKI signal models are "non-physical" and break down at high b-value, since the NTVD supports also tensors which are not positive definite, and the DKI signal decay is not the CF of a probability distribution.
3 Sufficiency of rank-1 and 2 b-matrices for diffusion tensor covariance estimation is spanned by tensors of the form where ⊗ denotes outer product, and by symmetrizing it follows that the tensors u ⊗ u ⊗ u ⊗ u with u ∈ R d span the subspace of totally symmetric 4th-order tensors.
Proof. By the spectral decomposition, every T ∈ T is a linear combination of at most (d+1)d 2 4th-order tensors of the form It can be observed that Equation (28) is a sum of outer squares of rank-1 and rank-2 b-matrices.
Thus, rank-3 b-matrices are not necessary and rank-1 b-matrices are not enough to estimate the covariance of the diffusion tensor in a multiple pulsed gradient experiment.

Well-posedness of signal inversion
The CNTVD is the probability distribution with given mean, D , and covariance, D, D , which is maximizing the entropy (relative to the volume measure) and therefore it is determined by its first and second moments. It is also shown in Theorem 1 that the mean and covariance determine uniquely the parameters µ, Σ of the CNTVD.
In Section 3 of the supplementary material we showed diffusion tensor covariance as fourth order partially symmetric tensors spanned by tensors, b ⊗ b, where b matrices are of rank 1 and 2. Below we show that from acquisitions with b-matrices of rank 1 and 2 we can unambiguously reconstruct the first and second moments of any diffusion tensor distribution, CNTVD included.
Consider the log-signal, log S(b) = log exp −Trace(bD) with a rank-1 b-matrix, εu u, where u is an unit vector and take its derivative at ε = 0 using the cumulant expansion (Equation (20), Since (by spectral theorem) rank-1 matrices span the symmetric matrices, it follows that rank-1 b-matrices and zero b-value are enough to determine the mean tensor and S 0 .
Consider the log-signal, log S(εb), with scalar ε > 0, and non-negative matrix b, and take its second derivatives at ε = 0 using the cumulant expansion (Equation (20)), Now because in Section 3 Lemma 1 the outer products b ⊗ b with b of rank-1 and rank-2 span the space of fourth order covariance tensors, in principle data from rank-1 and rank-2 b-matrices are sufficient to estimate the covariance of DTD. Since the CNTVD is fully described by the mean and covariance tensors, the above proof shows that the signal inversion is well-posed using rank-1 and rank-2 b-matrices.

Model selection hierarchy
The diffusion weighted data is fitted with the following family of nested models for the mean [2] tensor and covariance tensor [3] in the decreasing order of model parsimony as shown below. The    The size heterogeneity measure (V size ) is evaluated by mixing "RandomShapes" and "Ran-domShapes Resized" class of diffusion tensors. The shape heterogeneity measure (V shape ) is evaluated by mixing Oblate and Sphere class of diffusion tensors. The orientation heterogeneity measure (V orient ) is evaluated by mixing "RandomShapes" and "RandomShapes Rotated" class of diffusion tensors. Finally, "Oblate" and "ProlateRand" class of diffusion tensors are mixed together to evaluate the efficacy of the stains in filtering size, shape and orientation heterogeneities in the resulting mixture. The mixing fraction is gradually increased from 0 (all tensors belonging to the first group) to 1 (all tensors belonging to the second group) for controlled mixing of two groups.
The results of the mixing is shown in Figure 1 below. It can be observed that in the first case ("RandomShapes" vs "RandomShapes Resized"), all the stains except the size heterogeneity stain change with mixing fraction. The trend observed in the size heterogeneity stain was also parabolic with maximum attained at mixing fraction equal to 0.5 corresponding to the mixture with maximum entropy which reflects the ability of the stain to accurately measure the extent of size mixing of diffusion tensors. The shape heterogeneity stain was uniformly high due to the randomizing of shapes of both the mixing population. The shape heterogeneity measure in the second case ("Oblate" vs "Sphere") also exhibited a parabolic trend with mixing fraction but zero value for size and orientation heterogeneity measures as expected. The orientation heterogeneity measure in the third case ("RandomShapes" and "RandomShapes Rotated) also exhibited a parabolic trend with mixing fraction while the shape measure was uniform and above zero due to the presence of shape heterogeneity in both the mixing populations. The fourth case clearly showed parabolic trends for size and shape heterogeneity measures and approximately linear trend for the orientation measure as expected in the fully heterogeneous mixture. Oblate vs ProlateRand V orient V shape V size Figure 1: Efficacy of size (V size ), shape (V shape ) and orientation (V orient ) heterogeneity stains in independently capturing the desired heterogeneity. The plot of heterogeneity measures versus the mixing fraction computed for mixture of two groups of : (First column) randomly shaped tensors of same orientation but different sizes, (Second column) oblate and spherical tensors of same size, (Third column) randomly shaped tensors of same size but rotated by 90 • around the axis, [1,1,1], and (Fourth column) oblate and randomly rotated prolates of different sizes. The above graphics was generated using the matplotlib package [4]