Probing tissue microstructure by diffusion skewness tensor imaging

Probing the cellular structure of in vivo biological tissue is a fundamental problem in biomedical imaging and medical science. This work introduces an approach for analyzing diffusion magnetic resonance imaging data acquired by the novel tensor-valued encoding technique for characterizing tissue microstructure. Our approach first uses a signal model to estimate the variance and skewness of the distribution of apparent diffusion tensors modeling the underlying tissue. Then several novel imaging indices, such as weighted microscopic anisotropy and microscopic skewness, are derived to characterize different ensembles of diffusion processes that are indistinguishable by existing techniques. The contributions of this work also include a theoretical proof that shows that, to estimate the skewness of a diffusion tensor distribution, the encoding protocol needs to include full-rank tensor diffusion encoding. This proof provides a guideline for the application of this technique. The properties of the proposed indices are illustrated using both synthetic data and in vivo data acquired from a human brain.

Noninvasive probing of cellular structure of biological tissue in vivo is a fundamental problem in medical science. Modern medical imaging modalities, such as computed tomography (CT) and magnetic resonance imaging (MRI), provide imaging data with spatial resolution at the scale of millimeter (mm). But mental disorders and other brain diseases are related to alternations in the cellular microstructure at the micrometer ( µ m) scale without any gross effect at the macroscopic scale. However, no in vivo imaging technique is currently available to directly probe microscopic cellular arrangements in the human brain. Diffusion MRI (dMRI) is a modality to indirectly characterizes the microscopic cellular arrangements via the diffusion trajectories of water molecules 1,2 .
The diffusion coefficient of water at human body temperature is approximately 3 µm 2 /ms . In biological tissue, e.g. human brains, the displacement of water molecules is restricted or hindered by cellular membranes. The measured diffusion coefficients using dMRI, which is usually called the apparent diffusivity, depend on the microstructure of tissue as well as the experimental parameters such as the diffusion time. Within the typical 10 to 100 ms diffusion time used in dMRI, the radius of the diffusion trajectory of a water molecule is usually in the order of several µ m, similar to the size of cells. Because of microstructural heterogeneity, water molecules from different tissue components within one voxel in dMRI have different apparent diffusivity. Moreover, the diffusion trajectories can be characterized by different principal diffusion directions and different degrees of anisotropy. Thus, characterizing the statistical properties of all diffusion trajectories provide information to indirectly assess the property of tissue microstructure. Furthermore, dMRI can be performed with various diffusion encoding waveforms to sensitize the MR signal to different properties of diffusion trajectories. Then, tissue microstructure can be estimated by using a suitable analysis method for dMRI data. For example, the apparent diffusion tensor is a standard technique to characterize orientation-dependent diffusion which is related to the underlying cellular or axonal directions 3 . But this simplistic model only reflects the diffusion tensor of the ensemble average process with no information provided about the underlying variance and other high-order moments of the distribution of diffusivity of all water molecules that are useful to characterize the tissue heterogeneity.
In fact, this limitation is not only caused by the signal model but also the acquisition protocol. Specifically, dMRI data acquired by the standard linear-encoding sequence is not sufficient to characterize the variance of the distribution of apparent diffusion tensors 4 . To overcome this limitation, more advanced diffusion encoding waveforms, such as double-diffusion encoding 5,6 , q-MAS isotropic diffusion encoding 7 , and more general q-trajectory encoding (QTE) 4 www.nature.com/scientificreports/ tensor distributions, the relationship between the diffusion encoding sequences and higher-order moments, e.g. skewness, remains unclear. Currently, no method is available to use the high-order moments to analyze the heterogeneity of apparent diffusion tensors within one voxel. This work introduces an approach for modeling and analyzing dMRI data acquired using novel QTE waveforms in order to derive several novel imaging indices based on the skewness of the diffusion tensor distribution. These indices are able to distinguish different distributions of diffusion tensors that cannot be set apart by current approaches, providing new imaging measures of tissue microstructure. Moreover, this work also presents a theoretical proof that full-rank b-tensors are necessary to estimate the skewness tensor. The feasibility of the proposed approach is illustrated using an in vivo dataset acquired from a human brain.  Fig. 1 all have the same mean isotropic diffusion tensor. As a result, the three DTDs cannot be distinguished by any indices derived from the mean tensor, including the fractional anisotropy (FA) and the skewness (SK) see Methods. The microscopic anisotropy, µFA , has been proposed to characterize more specific information of the DTDs using the underlying variances 4,7 , see Methods. But the µFA of the three DTDs shown in Fig. 1 have similar µFA values, indicating the limited sensitivity and specificity of µFA for characterizing tissue microstructure.

Experiment and results
By using the third-order moment of the DTDs, we derived the microscopic skewness µSK and the weighted microscopic anisotropy, i.e. µFA fast and µFA slow , to characterize the statistical property of the DTDs, see Methods. In particular, µSK is sensitive to the shape of the microscopic diffusion tensors. In Fig. 1, DTD 1 and DTD 2 have negative and positive µSK values, respectively, which reflect the underlying oblate and prolate ellipsoids for tensor representations. On the other hand, µFA fast and µFA slow are sensitive to the microscopic anisotropy of diffusion tensors with fast and slow diffusivity. DTD 3 includes isotropic tensors with relatively fast diffusivity and anisotropic tensors with relatively slow diffusivity. As a result, the underlying µFA slow is higher than µFA fast . But µFA fast and µFA slow in DTD 1 and DTD 2 are not different, since the underlying tensors are of the same shape and size.
Diffusion encoding waveforms. The proposed statistical indices have been applied to analyze an in vivo dataset of acquired from a human brain using a clinical 3T Siemens MAGNETOM Prisma scanner, see Methods for experimental parameters. The experiment was approved by a research ethics committee of Brigham and Women's Hospital, Boston, MA. The first and second rows of Fig. 2 illustrate the diffusion-encoding gradient trajectories corresponding to different shapes of b-tensors, where the waveforms along x, y and z axes are coded by different colors. In the third row of Fig. 2, the rotation directions of the waveforms are shown by the dots on the spheres, and the b-values, i.e. the trace of b-tensors, are indicated by the radius of the spheres. A total number of 513 b-tensors are used in this experiment to ensure that the proposed statistical indices are uniquely determined by the dMRI signals, see Methods.
Microstructural measures of a human brain. Figure 3 shows the estimated dMRI measures, including FA, µFA, µFA slow , µFA fast , SK, µSK and the corresponding standard deviation (STD) in an axial slice of the brain. The STD values were estimated using the residual bootstrap method proposed in 10 . The first row shows  Fig. 3d shows that µFA is more reliable in white-matter. The second row shows that µFA slow is higher than µFA fast especially in gray matter. It was shown in 11 that extra-axonal water may have slower diffusivity than the intra-axonal water. Thus, extra-axonal water may be also related to higher µFA . Figure 3f,h show that µFA slow has lower STD values than µFA mfast especially in white matter. The third row shows the SK , the micro-SK, µSK and their STD maps. The white-matter region in the yellow circle has more negative SK values than the region corresponding to the red circle, though the FA in the yellow circle is relative higher. Thus, the diffusion tensors in the yellow-circle region are more prolate and anisotropic than those from the redcircle regions. The difference in FA and SK between the two regions may be related to the number of crossingfibers 12,13 , the fiber dispersion 14 and/or the angle between crossing fibers 15 , which needs to be further examined. The STD of SK in Fig. 3j is relatively high in white-matter, especially in regions with crossing fibers. The very high µSK values in the boundary between white/gray matter and the cerebrospinal fluid (CSF) in Fig. 3k are related to high STD values shown in Fig. 3l. It indicates that the µSK in these regions are sensitive to measurement noise, which may be related to the underlying partial volume effects, i.e. both CSF and white/gray matter are contained in the underlying voxel. Parameter uncertainty is inflated in a band toward the posterior of the brain as shown in Fig. 3b,f,h,j. This is likely due to poor fat saturation causing a mixture of fat and water signals at this location, as previously described by 9 .

Accuracy and precision of diffusion indices.
To examine the reliability of the estimated indices, the 513 b-tensors were applied to simulate dMRI signals for the three DTDs shown in Fig. 1 with additive Gaussian noise. The dMRI signals were computed by sampling s(B) = �e −B:D � ρ + ν with B representing one of the 513 b-tensors and ν representing the zero-mean Gaussian noise. For a fixed standard deviation of Gaussian noise, the 513 dimensional dMRI signals were sampled 5000 times. Then, the same estimation methods as applied in the in vivo experiment were applied to estimate the µFA fast , µFA slow and µSK measures. To further examine the accuracy and precision of the indices for different DTD functions, the experimental procedure was applied to the DTD function shown in Fig. 1c with varying volume fraction of anisotropic components and the signal to noise ratio (SNR), i.e. the ratio between the mean signal and the standard deviation of noise, being fixed at SNR=30 similar to the SNR values of the in vivo data.  Figure 5 shows the estimated indices corresponding to the DTD shown in Fig. 1c with different volume fractions for the anisotropic components. Similar to Fig. 4, the solid red and blue lines show the median of the estimated values and the dashed lines show the underlying true values. The upper and lower boundaries of the shaded area correspond to 25% and 75% percentiles of the estimated measures. The true values stay within the shaded areas when the volume fraction of anisotropic components is higher than 0.4. But marked errors appear as the volume fraction decreases, indicating that partial volume effect may lead to biases in the diffusion indices. Moreover, the µSK has very high uncertainties indicating that it is challenging to be correctly estimated, which is consistent to the results of the gray matter region in Fig. 3l.   Figure 3. Illustration of dMRI indices and the corresponding uncertainty of an axial slice of dMRI volume from a human brain. The first row shows the FA , µFA and the corresponding standard deviation (STD). The second row illustrates the µFA slow and µFA fast measures, where µFA slow has relative higher values than µFA fast . The third row shows SK and µSK . SK has negative values in the white-matter region highlighted by the red and yellow circle whereas the corresponding µSK are positive, indicating that the underlying microscopic diffusion tensors have prolate shapes that are related to crossing fibers or fiber dispersion. Parameter uncertainty is inflated in a band toward the posterior of the brain. This is likely due to poor fat saturation causing a mixture of fat and water signals at this location, as previously described by 9 .

Discussion and conclusions
This work introduces a method for estimating the variance and skewness of a distribution of diffusion tensors as an approach for characterizing the microstructure of materials or biological tissue. The method is based on expanding the diffusion tensor probability distribution in its cumulant moments, and relies on the diffusion weighting with full-rank b-tensors. The skewness of diffusion tensor 16 , diffusion skewness 17 and high-order tensors [18][19][20] have been investigated for tissue microstructure estimation based on standard (rank-one) b-tensor diffusion encoding (SDE). The added benefit of the proposed method with more general diffusion encoding is a novel capacity to distinguish diffusion tensors based on their skewness on the microscopic level, disentangling prolate and oblate shapes, as well as separately quantifying the microscopic anisotropy of diffusion tensors with high and low isotropic diffusivity. Taken together these features can add additional information on tissue microstructure, expanding on lower order methodology like DKI 21 and QTI 4 , which have been used to characterize prostate cancer 22 and intracranial tumors 23 . We note that the proposed model ignores the intra-compartmental  www.nature.com/scientificreports/ non-Gaussian diffusion 24,25 and time-dependent diffusivity 26,27 , which can be further integrated to the model to improve the accuracy for microstructure estimation. We do not expect these effects are significant in the healthy brain especially in white matter 28,29 . Moreover, the proposed indices based on cumulant moments, similar to DKI and QTI, do not directly provide specific information about tissue microstructure such as intra-axonal/cellular volume fractions, though they are sensitive to microstructural changes. From studies using cumulant expansions up to the second order we know that there exists a tradeoff between trueness and precision, where the use of higher b-values typically increase the precision while being detrimental to the trueness 30 . Nevertheless, metrics provided by such approaches can be useful in radiology for example to improve glioma grading 31 . The trueness of the parameters proposed here certainly also depends on the specifics of the imaging protocol. Further studies are needed to carefully balance the trueness and precision in this context. We emphasize the need for further investigations of tissues of pathological conditions where this may provide relevant information. We speculate that it may resolve the equality of prolate and oblate solutions and being able to robustly detect the presence of planar structures, such as those conjectured to exist in white matter sheets 15 . Currently, a proof-of-concept is illustrated in a healthy human brain in vivo. Indices derived from the model will be validated using liquid crystal phantoms 32 and biological tissue in future works. In summary, this work has introduced an approach for modeling and analyzing diffusion MRI data acquired by the novel sampling scheme with a wide range of b-tensor shapes. The main contributions include several statistical indices based on the skewness tensor of DTD functions and a mathematical proof that full-rank b-tensors are needed to uniquely determine the skewness tensor.

Modeling of diffusion MRI signals.
In a diffusion experiment, on the application of a time-varying magnetic gradient field, g(t) ∈ R 3 for t ∈ [0, τ ] , the phase change of a rotating spin due to diffusive motion of water molecule during this time window is given by 33 : where γ is the gyromagnetic ratio, r x (t) denotes the trajectory of the displacement of a particle starting from x , v x (t) denotes the corresponding velocity process, which is formally defined as v x (t) =ṙ x (t) , q(t) = γ t 0 g(s)ds is the q-trajectory and T denotes the transpose. If the diffusion trajectory is assumed to be a three-dimensional Wiener process with constant diffusivity D during the period of diffusion time, then φ x (τ ) follows a zero-mean Gaussian distribution with the covariance equal to where B = τ 0 q(t) ⊗2 dt is the b-tensor. Then, the expected value of the corresponding MR signal is equal to E (e iφ x (T) ) = e −B:D .
The diffusion-weighted MR signal after being normalized by the non-diffusion-weighted signal, i.e. the baseline, is equal to the ensemble average of signals from all water molecules given by: where �·� p denotes the ensemble average with respect to the phase distribution function p(φ) . At long diffusion time scale when the travel distances of water molecules much longer than characteristic cellular sizes, the molecular diffusion processes can be approximated by Gaussian processes with location dependent diffusivity 34 . At long diffusion time, if the mixture of apparent diffusion tensors from different cellular compartments within a voxel follows the diffusion tensor distribution (DTD) function ρ(D) , then the diffusion-weighted signal is equal to 4,35 : A similar signal model based on scalar-valued diffusivity was investigated in 36 .
The cumulant expansion. The standard dMRI approaches probe the diffusion process using a linear diffusion gradient waveform which corresponds to a b-tensor of rank one. In this case, estimating the underlying DTD ρ(D) using an inverse Laplace transform is a highly ill-posed problem. However, advanced QTE sequences are able to probe ρ(D) using any b-tensor up to rank 3, providing information that cannot be measured by SDE 5 . The ensemble average signal �e −B:D � ρ can be approximated using the following cumulant expansion: where D ρ , C, S denote the first three cumulants and are given by: C = �(D − �D� ρ ) ⊗2 � ρ and S = �(D − �D� ρ ) ⊗3 � ρ . The notation ⊗ denotes the tensor product and " : ′′ defines the standard inner product between matrices. For example, B ⊗3 : S = [B ⊗3 ] ijk [S ] ijk using Einstein summation convention. The tensors D ρ , C and S have 6, 21 and 56 independent variables, respectively. These 83 variables can be estimated by solving the following linear system of equations:  Similarly, − → C and − → S are 21-and 56-dimensional column-vector representations for C and S , respectively. Clearly, the number of measurements, m, needs to be at least 83 to ensure that the measurement matrix ( m × 83 ) is non-singular. The 513 b-tensors used in the experiment provides a full-rank measurement matrix of size 513 × 83 . Next we show that non-singular b-tensors are needed to ensure that the measurement matrix is of full rank.
On full-rank b-tensors. In linear encoding sequences, the gradient sequences have varying magnitudes along a fixed direction. In this case, the corresponding b-tensors have rank one. But the corresponding 21-dimensional vector −→ B ⊗2 only spans a 15-dimensional subspace. Consequently, the variable − → C in the linear system (3) cannot be uniquely determined using standard SDE measurements 4 . However, the vectors −→ B ⊗2 corresponding to rank-2 b-tensors given by planar encoding sequences 5,6,37,38 or rank-3 tensors general QTE sequences 9 are able to generate full-rank 21-dimensional space, providing a unique solution to C . But b-tensors of rank two still have zero determinant, i.e. which implies that the corresponding elements in −→ B ⊗3 are linearly dependent. Therefore using only b-tensors of rank one or rank two is not sufficient to generate full-rank measurement matrices in Eq.
(3). This shows that a full-rank b-tensor is needed to ensure a unique solution to Eq. (3). The importance of rank-2 b-tensors in estimating the covariance of DTDs was considered in 5,8 . The above analysis further extends the theory to show the importance of full-rank b-tensors in estimating the skewness of DTDs.
Microscopic anisotropy and skewness. The estimated moments (or cumulants) can be used to derive scalar indices to characterize the underlying DTDs. In particular, the standard diffusion tensor imaging (DTI) technique uses the mean diffusion tensor D ρ to compute several classical indices including the mean diffusivity, MD = 1 3 trace(�D� ρ ) and the fractional anisotropy FA 3 . But these measures provide limited information about the underlying DTD. For example, the mean diffusion tensor of the three DTDs illustrated in Fig. 1 is proportional to the identity matrix with FA = 0, though the underlying DTD functions are different. One can look at the covariance tensor C , and compute the micro-anisotropy µFA of the DTDs given by: 3 I 6×6 with I 6×6 being the identify matrix of size 6 × 6 and E shear is equal to 4 However, in many scenarios, even the variance tensor is not enough to distinguish between the DTDs. For example, the three DTD functions in Fig. 1 have the same µFA measure though the underlying microscopic diffusion tensors are different. In this case, the third-order moments provide useful information to distinguish the underlying DTD functions. To this end, we propose two sets of indices, µFA fast , µFA slow and µSK , to characterize the DTDs using the third-order moments.
To introduce the first approach, we introduce a family of DTD functions defined by where f(D) is a given positive scalar-valued function of the diffusion tensor D. This function is chosen to scale the probability density so that the statistical property of the filtered DTD function ρ f (D) provides more specific information about the underlying distribution. The moments of the filtered DTD function ρ f (D) can be computed using: These equations can be used to estimate specific dMRI measures using the estimated moments for ρ fast (D) and ρ slow (D) . For the DTD functions shown in Fig. 1, the µFA measures for ρ fast (D) and ρ slow (D) , denoted by µFA fast and µFA slow , corresponding to DTD 1 and DTD 2 have the same values since the underlying diffusion tensors have the same trace. For DTD 3 , µFA fast is much lower than µFA slow since the underlying fast diffusion tensors are isotropic whereas the slow-diffusion components are anisotropic.
To introduce the microscopic skewness, we first extend the skewness measure for probability distribution functions to define the following measure of asymmetry of the distribution of eigenvalues of the mean diffusion tensor D : where D is assumed to be anisotropic so that the denominator is assumed nonzero. In above, E shear is a fully symmetric three-dimensional tensor, i.e. the value of [E shear ] ijk does not change by permuting the order of i, j, k, of size 6 × 6 × 6 with the following non-zero entries where I 4 = {1, 2}, I 5 = {1, 3} , I 6 = {2, 3} . The skewness measure SK is able to distinguish diffusion tensors that cannot be separated by FA values. For example, the diffusion tensors represented by oblate tensors have a negative SK while prolate tensors have positive SK . Following the definition of µFA , we define the following microscopic skewness measure: (7) f fast (D) = trace(D),   www.nature.com/scientificreports/ which characterizes the microscopic skewness of the diffusion tensors. To enhance robustness with respect to measurement noise, the �D ⊗2 � ρ : E shear term in the denominator on the right hand side of the above equation is replaced by �D ⊗2 � ρ : E shear + ǫ with ǫ = 0.03 µm 4 /ms 2 , which does not change the sign of µSK(ρ) . In particular, the µSK(ρ) for DTD 1 in Fig. 1 is negative while for DTD 2 and DTD 3 it is positive.
Experimental parameters of the in vivo dataset. The "free waveform" prototype pulse sequence 9 based on a diffusion-weighted spin-echo sequence with EPI readout was used to acquire MRI data from an adult male subject. Imaging parameters were: TE = 100 ms, TR = 2800 ms, field of view = 220 × 220 × 75 mm 3 , voxel size = 2.2 × 2.2 × 5 mm 3 , partial Fourier factor = 6/8 and in-plane acceleration factor = 2. The acquired images were corrected for motion and eddy current distortion using an extrapolation-based approach capable of correcting diffusion weighted volumes acquired with high b-values 39,40 . A total number of 513 image volumes were acquired using different QTE trajectories tailored to the MRI system by numerical optimization 41 with compensation for concomitant gradient effects 42 . The QTE waveforms were determined by two parameters b and b η , which characterize the shape of b-tensors 9,43 . The orientation of b-tensors was determined via rotation of the waveforms and b-values were changed by varying the magnitude of waveforms 44 . The diffusion encoding waveforms in the first two columns provide spherical and oblate b-tensors with full-rank, which are required to estimate the skewness of the underlying DTD. The 513 b-tensors in this experiment provide a full-rank measurement matrix so that the proposed statistical indices are uniquely determined by the dMRI signals.