Evaluation of non‐Gaussian diffusion in cardiac MRI

Purpose The diffusion tensor model assumes Gaussian diffusion and is widely applied in cardiac diffusion MRI. However, diffusion in biological tissue deviates from a Gaussian profile as a result of hindrance and restriction from cell and tissue microstructure, and may be quantified better by non‐Gaussian modeling. The aim of this study was to investigate non‐Gaussian diffusion in healthy and hypertrophic hearts. Methods Thirteen rat hearts (five healthy, four sham, four hypertrophic) were imaged ex vivo. Diffusion‐weighted images were acquired at b‐values up to 10,000 s/mm2. Models of diffusion were fit to the data and ranked based on the Akaike information criterion. Results The diffusion tensor was ranked best at b‐values up to 2000 s/mm2 but reflected the signal poorly in the high b‐value regime, in which the best model was a non‐Gaussian “beta distribution” model. Although there was considerable overlap in apparent diffusivities between the healthy, sham, and hypertrophic hearts, diffusion kurtosis and skewness in the hypertrophic hearts were more than 20% higher in the sheetlet and sheetlet‐normal directions. Conclusion Non‐Gaussian diffusion models have a higher sensitivity for the detection of hypertrophy compared with the Gaussian model. In particular, diffusion kurtosis may serve as a useful biomarker for characterization of disease and remodeling in the heart. Magn Reson Med 78:1174–1186, 2017. © 2016 International Society for Magnetic Resonance in Medicine.


INTRODUCTION
Diffusion MRI is used to noninvasively provide information about tissue microstructure (1). Diffusion tensor imaging (DTI) is a widely used method in which a minimum of six diffusion-weighted (DW) images, in addition to a reference image without diffusion weighting, are used to derive a diffusion tensor (2). In cardiac MRI, the primary, secondary, and tertiary eigenvectors of the diffusion tensor have been shown to correspond to the socalled fiber, sheetlet, and sheetlet-normal orientations within the myocardium, respectively (3,4). Indices derived from DTI such as apparent diffusion coefficient (ADC) and fractional anisotropy (FA) can be linked broadly to a range of tissue properties, including water molecule mobility and coherence of regionally prevailing cell orientation (fiber orientation). In the heart, these indices have been used to assess remodeling of tissue structure in disease states such as hypertrophic cardiomyopathy (5) and after myocardial infarction (6). The diffusion tensor model assumes that the displacement profile of water diffusion follows a Gaussian distribution. However, it is well established that diffusion in biological tissue deviates from a Gaussian profile as a result of hindrance and restriction from cell and tissue microstructure (7). To address this, a number of models have been proposed to more accurately describe the non-Gaussian behavior of diffusion in tissue, including the diffusion kurtosis model (8), biexponential models (9), stretched exponential models (10), and statistical models (11)(12)(13). These models have been widely applied in brain MRI, and it has been shown that models of diffusion kurtosis can act as a biomarker for clinical development of Alzheimer's disease (14), provide improved characterization of microstructure in the developing brain (15), and improve characterization of tumors (16) compared with diffusion tensor modeling. Compartmental models of diffusion have also been successful in estimating axon diameter, density, and dispersion (17).
To date, non-Gaussian analysis of cardiac MRI data has been limited. A bi-exponential model has been investigated in perfused rat, rabbit, and guinea pig hearts, with the fast and slow components associated with contributions from perfusion and diffusion, respectively (3,18,19). It has also been shown that non-monoexponential diffusion is present in both healthy and infarcted fixed rabbit hearts (20). However, this study was limited to quantifying deviations from the diffusion tensor and did not use models of non-Gaussian diffusion.
Non-Gaussian diffusion models may offer metrics that are more sensitive to the presence of restrictions such as cell membranes and organelles. This could be particularly relevant in disease states such as hypertrophy, in which proliferation and enlargement of cardiomyocytes leads to changes in the distribution of cellular restrictions to water diffusion (21). Transverse aortic constriction (TAC) is widely used to surgically generate animal models of heart failure (22). Partial constriction of the transverse aorta via a metal clip or suture results in a rapid increase in left ventricular load with consequent development of hypertrophy, ultimately leading to heart failure.
The aim of this study was to systematically investigate non-Gaussian diffusion in healthy and hypertrophic fixed rat hearts. A conventional diffusion tensor model as well as a number of non-Gaussian models were fit to DW images of rat hearts, and ranked based on the Akaike information criterion (AIC) (23). Next, the non-Gaussian diffusion of healthy hearts was compared with hypertrophic hearts. We hypothesize that non-Gaussian diffusion modeling provides more sensitive biomarkers in hypertrophic hearts than diffusion tensor modeling. A preliminary version of this work was presented previously (24).

THEORY
In a standard pulsed-field-gradient sequence with rectangular gradient pulses of strength G (in T/m), diffusion gradient duration d (in s), and diffusion time D (in s), the diffusion weighting is quantified by the b-value, given where c is the gyromagnetic ratio (in rad/s/T). Models of diffusion are used to parameterize the relationship between the b-value and the measured signal attenuation. Assuming a voxel contains a spectrum of diffusion environments, each exhibiting Gaussian diffusion, the measured signal is given by where PðDÞ describes the probability density function (PDF) of diffusivity D. Various distributions have been proposed for PðDÞ, including the truncated Gaussian (13), gamma (11,25,26), log-normal (12,27), and beta (26) distributions. For a given PðDÞ, the mean diffusivity is given by the expectation of D, E½D. Excess kurtosis (referred to hereafter as kurtosis) is given by the normalized variance of the diffusivity profile as follows: The displacement profile is obtained via the inverse Fourier transform of the signal with respect to the wave vector q ¼ gdg 2p (28). Figure 1 presents the signal attenuation, diffusivity profile, and displacement profile of selected models. These models are primarily sourced from the review article of Yablonskiy and Sukstanskii (26) and are described in greater detail hereafter. Additional details about mathematical functions required for the models are given in the appendix.

Monoexponential Model
The well-known monoexponential model makes the assumption that a voxel contains a single compartment exhibiting unrestricted diffusion. It is derived from the Taylor series of the relationship between the b-value and logarithm of the MRI signal SðbÞ as follows: ln SðbÞ Sð0Þ ¼ ÀbD app þ Oðb 2 Þ: [3] Truncating this series at the first term yields the monoexponential model Evaluation of Non-Gaussian Diffusion in Cardiac MRI SðbÞ Sð0Þ ¼ expðÀbD app Þ; [4] where D app is the apparent diffusion coefficient (in mm 2 / s). The PDF associated with the monoexponential model is a delta function at D app .

Stretched Exponential Model
The stretched exponential model was first introduced in 1854 to describe the discharge of a capacitor (29). It was proposed in the context of diffusion MRI by Bennett et al. (10) and later linked to the concept of fractal dimension (30). This model is defined as where a is a dimensionless stretching index, 0 < a 1 and D s is the "stretch-adjusted" diffusivity. Neither the mean diffusivity nor the diffusion kurtosis can be derived from this model because b a cannot be described by a Taylor series for noninteger values of a. However, a decreased stretching index a may intuitively be interpreted as having greater deviation from a Gaussian displacement profile, which is associated with higher kurtosis. A decreased stretching index has also been associated with more crowded microstructure and more complex water excursions (30). The diffusivity PDF does not have an analytical form but can be computed numerically (31).

Diffusion Kurtosis Model
The diffusion kurtosis (DK) model (8) is derived from the Taylor series expansion of ln SðbÞ Sð0Þ truncated at the second-order term and is defined as follows: The mean diffusivity and kurtosis of this model are given by D DK and K DK , respectively. The PDF of the DK model is a Gaussian distribution with mean D DK and variance 1 3 D 2 DK K DK . Thus, kurtosis is a measure of the heterogeneity of the diffusion environment, where higher values of kurtosis indicate a wider spread of apparent diffusivities within a voxel. The decay curve is a quadratic on a log-scale, with a minimum at b ¼ 3 DDKKDK . The selection of b-values is important in diffusion kurtosis imaging; the maximum b-value must be high enough that that the Oðb 2 Þ is measurable in the presence of noise, but low enough that the truncated Oðb 3 Þ term is negligible.

Biexponential Model
The biexponential model (9) assumes that the diffusion signal attenuation can be attributed to two populations, described as a fast and slow component, each of which exhibits Gaussian diffusion.

SðbÞ Sð0Þ
¼ v expðÀD fast bÞ þ ð1 À vÞexpðÀD slow bÞ: [7] Parameter v is the volume fraction of the fast component. Mean diffusivity and kurtosis may be computed from the biexponential model as follows: The biexponential model corresponds to a PDF consisting of a weighted sum of two delta functions.

Truncated Gaussian Distribution Model
The truncated Gaussian distribution (13) is appealing because it approximates the DK model at low b-values but is monotonically decreasing at arbitrarily high bvalues. The PDF is defined as where D m and s (both in mm 2 /s) control the location and scale of the distribution, respectively, and A is a normalizing coefficient. The normalizing coefficient is dependent on both D m and s. The observed signal from this model is given by where erfðÁÞ denotes the error function (see Eq. [A1]).
Mean diffusivity, D, is given by D m þ s ffiffi ffi , and kurtosis is given by 3

Gamma Distribution Model
The gamma distribution (11,25,26) allows more variation in the shape, and in particular the skewness, than the truncated Gaussian model. The PDF is defined as where G denotes the gamma function (see Eq. [A2]). Parameters k and u control the shape of the distribution. The observed signal is given by Mean diffusivity is given by ku, and kurtosis is given by 3 k . The skewness of this distribution is given by 2 ffiffi ffi k p .

Beta Distribution Model
A limitation of the truncated Gaussian and gamma distributions is the presence of positive tails with diffusivity greater than that of unimpeded water. The beta distribution model (26) addresses this by imposing an upper limit D max on the PDF as follows: ; D2ð0; D max Þ ; [14] where B denotes the beta function (see Eq. [A3]) and parameters a and b control the shape of the distribution. The observed signal is given by where M is the confluent hypergeometric function (see Eq. [A4]). Because D max is generally not known, it can be treated as another parameter. Mean diffusivity and kurtosis are given by a aþb D max and 3 b a ðaþbþ1Þ , respectively. Skewness is given by 2ðbÀaÞ ffiffiffiffiffiffiffiffiffiffiffi ffi

Tissue Preparation
Experimental , and stored at 4 C. Prior to imaging, the hearts were rinsed three times in phosphatebuffered saline with 2 mM gadolinium chelate and embedded in 1% agarose gel (Web Scientific, Crewe, UK) made with phosphate-buffered saline containing 2 mM gadolinium chelate. Gel was used to retain sample geometric stability, immobility, and hydration, and gadolinium chelate shortened T 1 and increased signal-tonoise ratio (SNR) efficiency. Eight additional male Sprague-Dawley rats (body weight, 206 6 5 g) were divided into sham (n ¼ 4) and TAC (n ¼ 4) groups. Animal studies for the sham and TAC groups followed National Institutes Health guidelines and were approved by the University of California San Diego Institutional Animal Care and Use Committee. Animals were anesthetized to an appropriate level with 1.25%-2% inhaled isoflurane, which was confirmed via toe pinch. A lateral thoracotomy exposed the aortic arch. Aortic stenosis was induced in animals in the TAC group by means of a ligature hemo-clip on the transverse aorta between the right innominate and left carotid arteries with a final, constricted width of 0.5 mm. Sham controls underwent thoracotomy without application of the clip. The surgical site incisions were closed in layers, and the animal was allowed to recover. Four weeks after surgery, M-mode echocardiographic measurements were made, left ventricular pressures were recorded, and hearts were excised, arrested, fixed, and embedded for imaging in the same manner as the healthy hearts. Imaging was performed approximately 2 wk after fixation.

Imaging
Diffusion spectrum imaging data were acquired with qspace sampled on a three-dimensional (3D) Cartesian grid. A nonselective 3D fast spin echo sequence was used. Imaging was performed on a 9.4T horizontal bore MRI scanner (Agilent, Santa Clara, California, USA) with shielded gradients (max gradient strength ¼ 1 T/m, rise time ¼ 130 ms), and a transmit/receive quadrature-driven birdcage coil (inner diameter ¼ 20 mm; Rapid Biomedical, Rimpar, Germany).
For the healthy hearts, the following acquisition parameters were used: pulse repetition time ¼ 250 ms; echo time ¼ 15 ms; echo spacing ¼ 4 ms; echo train length ¼ 8; matrix ¼ 100 Â 80 Â 80; field of view ¼ 20 Â 16 Â 16 mm; isotropic resolution ¼ 200 mm; number of non-DW images ¼ 4; number of DW directions ¼ 257; b max ¼ 10,000 s/mm 2 ; diffusion duration d ¼ 5 ms; diffusion time D ¼ 9 ms; acquisition time ¼ 14 h, 30 min. The images were acquired at 20.2 C. Systematic temperature bias was minimized by pseudo-randomizing the DW acquisition such that consecutive scans were acquired with widely different b-values (32). Diffusion MRI experiments were followed by high-resolution anatomical imaging using the following parameters: pulse repetition time ¼ 20 ms; echo time ¼ 4 ms; flip angle ¼ 30 ; matrix ¼ 600 Â 480 Â 480; field of view ¼ 20 Â 16 Â 16 mm; isotropic resolution ¼ 33 mm; acquisition time ¼ 10 h, 14 min. The acquisition protocol for the sham and TAC hearts was identical to that for the healthy hearts, with the following exceptions: matrix ¼ 120 Â 80 Â 80; field of view ¼ 21.6 Â 14.4 Â 14.4 mm; isotropic resolution ¼ 180 mm; number of DW directions ¼ 514; acquisition time ¼ 28:46 h. Symmetric q-space acquisition and higher imaging resolution was used in the latter protocol due to greater availability of scan time.

Data Analysis
The noise distribution of the magnitude of MR signals with complex Gaussian noise is Rician. The non-zero mean of the noise results in biased estimates of parameters when fitting using least-squares minimization. This is particularly apparent in data with low SNR, as is common in high b-value imaging. To mitigate this bias, the phase was removed, yielding a zero-mean Gaussian noise distribution (33). Specifically, the mean of the complex non-DW images was computed, and the phase was subtracted from both the non-DW and DW images. The real component of the resulting images was extracted. These images can be assumed to have normally distributed noise and therefore can accommodate model fitting using least squares minimization algorithms. The noise level was estimated using the method described previously (34).
The diffusion tensor model was fit to the phasecorrected data. The eigenvalues of tensor D are given by k ¼ ½l 1 l 2 l 3 . The mean ADC and FA are defined as Helix, transverse, and sheetlet angle maps were computed from the eigenvectors of the diffusion tensor, using the definitions of angles by Hales et al. (35). The angles were defined relative to a local coordinate system as described by Teh et al. (32) to mitigate biases in angle maps arising from local tissue deformations. The helix angle is the angle subtended by the projection of the primary eigenvector onto the circumferential-longitudinal plane and the short-axis plane. The transverse angle is the angle subtended by the projection of the primary eigenvector onto the short-axis plane and the circumferential-longitudinal plane. The sheetlet angle is the angle subtended by the projection of the tertiary eigenvector onto the longitudinal-radial plane and the long axis. We follow the approach of De Santis et al. (36) in describing the observed attenuation from non-Gaussian models as the product of the attenuation in the direction of each of the three eigenvectors, as follows: where b i is the magnitude of the diffusion weighting b in the direction of the ith eigenvector of D and f is given by the right-hand side of Eqs. 4-7, 11, 13, and 15. This approach assumes that non-Gaussian diffusion shares the reference frame of the diffusion tensor model. Thus, each of the non-Gaussian models required six additional parameters to define the reference frame. A list of models in this study and total number of parameters for each is given in Table 1.
Each of the models was fit to all data sets (healthy, sham, and TAC) using nonlinear least squares regression. All models were fit in MATLAB R2013A (MathWorks, Natick, Massachusetts, USA) using a trust-regionreflective algorithm (37), with positivity constraints enforced where appropriate. The source code can be obtained through the gerardus project (https://github. com/vigente/gerardus/tree/papers). The initial values of the diffusion tensor and diffusion kurtosis model parameters were obtained using linear least squares regression on the logarithm of the signal. The input data to the DK model was restricted to a maximum b-value of 5000 s/ mm 2 (25).
Initial values of biexponential model parameters were determined by fitting a monoexponential model to low (<1000 s/mm 2 ) and high (>6000 s/mm 2 ) b-values to estimate the fast and slow diffusion rate, respectively (38). The stretched exponential model was initialized with the fit from the diffusion tensor model (i.e., a ¼ 1). In the case of the truncated Gaussian, gamma, and beta distribution models, initial estimates of parameters were derived using the mean diffusivity and kurtosis of the DK model. The parameter controlling the maximum diffusion coefficient in the beta distribution model, D max , was initialized with the diffusion coefficient of pure water at room temperature, 2.3 Â 10 À3 mm 2 /s (39).

Region of Interest Selection
A 3D region of interest (ROI) was manually traced on diffusion kurtosis maps in each heart, with reference to high-resolution anatomical images (Fig. 2a). In each case, the ROI was located in the myocardium of the left ventricle lateral wall. Care was taken to avoid regions containing buffer, gel, or residual blood, as the kurtosis in these regions was artificially increased, as shown in Figure 2b. The ROIs had a volume of at least 1.6 mm 3 (247 6 53 voxels).

Evaluation Using AIC c
Complex models having a larger number of fitting parameters provide lower residual errors, but may overfit the data, resulting in high parameter variance in the presence of noise. Conversely, simpler models with fewer parameters return lower parameter variance, but may fail to recover relevant information in the data.
This trade-off can be quantified using the AIC (23), which is based on information theory. The model selected by the AIC is the one that minimizes information loss. The AIC was selected for model comparison because it provides a quantitative estimate of the amount of information lost by different models and does not require the arbitrary selection of a cutoff. Performance was compared using the AIC with a correction for finite sample sizes (40), given by where e 2 is the mean squared error, n is the number of DW images, and P is the number of model parameters (including e 2 ). Thus, P was equal to 8 for the DTI model; 14 for the stretched exponential, diffusion kurtosis, truncated Gaussian, and gamma distribution models; and 15 for the biexponential and beta distribution models. Because AIC c is a relative measure, the difference between the AIC c of each model and the minimum AIC c was recorded. The relative likelihood of model i (i.e. the probability that this model minimises the estimated information loss) is given by where AIC min c refers to the minimum AIC c of the models considered.
To demonstrate b-value dependence, all models were fit to the data in the ROIs in each of the nine hearts, with diffusion data restricted to b-values lower than the maximum b-value. The maximum b-value was varied between 1200 s/mm 2 and 10,000 s/mm 2 in intervals of 400 s/mm 2 . AIC c and relative likelihood was computed in each case. The mean of the relative likelihood over all hearts in each of the three categories (healthy, sham, and TAC) is presented.

Pressure Overload
Aortic pressures measured in the TAC group before fixation verified the pressure overload, with the TAC animals having a greater left ventricular peak pressure compared with sham (198 6 43 mm Hg versus 120 6 10 mm Hg). This resulted in an increase in heart weight to body weight ratio (4.8 6 0.2 g/kg versus 3.5 6 0.2 g/kg). Furthermore, the interventricular septum thickness was higher in TAC versus sham (1. 46   Helix angles indicate the well-known continuous transition between left-handed and right-handed cell orientations. Transverse angles close to 0 indicate that the cells   FIG. 2. High-resolution anatomical image (a) and kurtosis in the tertiary eigenvector (b) in heart 6 (sham). High kurtosis correlates well with regions in the anatomical image that contain buffer. ROIs were traced in hypointense regions on kurtosis maps, with reference to the anatomical images. The ROIs were used for all comparisons between the healthy, sham, and TAC hearts.

Evaluation of Non-Gaussian Diffusion in Cardiac MRI
have a predominantly circumferential orientation when projected onto the short-axis plane. Sheetlet angle maps display distinct clusters of voxels with similar sheetlet angles.
Box-and-whisker plots of parameters derived from the diffusion tensor model in the ROIs are presented in Figure 4. On average, the diffusivity in the secondary and tertiary eigenvectors was lower in the TAC hearts compared with the sham hearts, resulting in increased fractional anisotropy. However, there is a large overlap in all DTI-derived parameters in the healthy, sham, and TAC hearts. The interquartile ranges of the normal hearts are greater than that of the sham and TAC hearts, most likely because of the smaller number of samples in the normal hearts (29 unique samples in q-space versus 41 in the sham and TAC hearts with b-values < 2000 s/ mm 2 ) leading to lower precision.

Model Selection Using AIC c
Model selection was highly dependent on the maximum applied b-value: Figure 5 presents the relationship between model fitting performance and the maximum b-value employed during fitting. In all hearts, the beta distribution model yielded the model with the lowest root-mean-square error (RMSE) and highest relative likelihood at maximum b-values up to 10,000 s/mm 2 (i.e., when all data were considered). At maximum b-values between approximately 3000 s/mm 2 and 7000 s/mm 2 , the truncated Gaussian model yields the highest relative likelihood in the healthy and TAC hearts, despite having a higher RMSE than the beta or biexponential models. In the sham hearts, the DK model outperformed the truncated Gaussian model for maximum b-values between 3800 s/mm 2 and 5400 s/mm 2 , although the difference between the two models is small. The stretched exponential model yielded the highest relative likelihood for maximum bvalues between approximately 3000 and 3800 s/mm 2 in the sham and TAC hearts. Finally, in all hearts, the conventional diffusion tensor model yielded the highest relative likelihood for maximum b-values up to approximately 2000 s/mm 2 . Figure 6 presents maps of the model with the highest relative likelihood at maximum b-values of 2000 and 10,000 s/mm 2 . At 2000 s/mm 2 , the diffusion tensor model yields the highest relative likelihood for the majority of voxels in gel, buffer, and the myocardium. The stretched exponential model yields the highest relative likelihood in most voxels on the interfaces between myocardium and gel/buffer. At 10,000 s/mm 2 , the beta model most commonly has the highest relative likelihood in the myocardium. The truncated Gaussian model has the highest relative likelihood in a large number of voxels, particularly in the right ventricle of the healthy heart. The biexponential model has the highest relative likelihood in some myocardial voxels, particularly in the sham heart, and very frequently has the highest relative likelihood on the interfaces between myocardium and gel/buffer. In the voxels on these interfaces, the fast component of the model corresponded to the gel/buffer with diffusivity of approximately 2.3 Â 10 À3 mm 2 /s. Figure 7 presents box-and-whisker plots of kurtosis and skewness generated from the beta distribution model in healthy, sham and TAC hearts, using all b-values up to 10,000 s/mm 2 . These data are also presented in Table 2. The beta model was selected due to its superior performance when fitting these data, as described above. Both The sham hearts had similar kurtosis and skewness values to the healthy hearts in all three eigenvectors. Two-tailed Mann-Whitney U tests were performed on the mean of model parameters in each of the ROIs in the sham and TAC hearts, at a 5% significance level. The kurtosis and skewness of the TAC hearts did not differ from that of the healthy or sham hearts in the primary eigenvector, but was significantly higher in the secondary and tertiary eigenvector directions. On average, kurtosis in TAC hearts was 31% higher than in sham hearts in the secondary eigenvector direction (P < 0.05), and 20% higher in the tertiary eigenvector direction (P < 0.05). Similarly, the skewness of the TAC hearts was 86% higher in the TAC hearts than in the sham hearts in the secondary eigenvector direction (P < 0.05), and 81% higher in the tertiary eigenvector direction (P < 0.05). Although the mean maximum diffusivity in the TAC hearts was slightly higher (1.50 mm 2 /s) than in the healthy (1.46 mm 2 /s) or sham (1.43 mm 2 /s) hearts, there was a large degree of overlap between the samples in the three categories.

Biexponential Modeling
Biexponential model parameters are presented in Table 2. The mean volume fraction of fast diffusing components was significantly lower in ROIs in the TAC hearts (0.82 6 0.01) than in the healthy (0.86 6 0.00) or sham (0.86 6 0.01) hearts (P < 0.05). On average, the fast diffusing components in the TAC hearts had higher diffusivity in v 1 and lower diffusivity in v 2 and v 3 than the healthy and sham hearts. However, the variance across the hearts in each category was relatively high. Similarly, differences between the TAC and healthy/sham hearts in the slow diffusing component were small relative to the interheart variance. Figure 8 presents the PDFs derived from the beta distribution model in the three eigenvectors. These were generated from averaged parameters in the ROIs of healthy, sham, and TAC hearts. The distributions in the primary eigenvector contain sharp peaks at 1.3-1.5 Â 10 À3 mm 2 /s. The distributions in the secondary and tertiary eigenvectors contain peaks at the same location, but with a broader spectrum (i.e., higher kurtosis) and greater contributions from low diffusivities (i.e., more positive skewness).

Diffusivity PDFs
The diffusivity profiles of TAC hearts in the secondary and tertiary eigenvectors are characterized by broader profiles than the healthy or sham hearts, indicating greater non-Gaussianity. This is in agreement with the higher volume fraction of the slow diffusion compartment in the biexponential model.

DISCUSSION
In this study, we investigated the non-Gaussian behavior of diffusion in ex vivo rat hearts using MRI. Hypertrophic hearts were compared with healthy and sham hearts on the hypothesis that non-Gaussian parameters FIG. 4. Box-and-whisker plots of parameters derived from the diffusion tensor model. Although the mean of the TAC hearts was lower than that of the sham hearts in l 2 and l 3 , there is considerable overlap, particularly between hearts 7 and 8. As a result, there is poor separation between hearts in terms of mean ADC and FA. will be more sensitive to hypertrophy than parameters derived from the Gaussian model.
In all hearts, greater non-Gaussianity was observed in the directions of the secondary and tertiary eigenvectors. This may arise from greater restrictions to water diffusion perpendicular to cardiomyocyte orientations. Non-Gaussianity was highest in the direction of the tertiary eigenvector, which may arise from increased heterogeneity in diffusion environments in the sheetlet-normal direction. In contrast, the profile of diffusivities in the primary eigenvector was narrower, indicating greater coherence of diffusion and a less restricted environment. These results are in agreement with a previous study on ex vivo pig hearts (41). Nguyen et al. (42) found higher ADC in putative fibrotic regions, as identified with ADC or extracellular volume (ECV) thresholds, in hypertrophic cardiomyopathy patients. This contrasts with our finding that ADC was similar in the sham and TAC groups. It is important to note here that Nguyen et al. focused on the clinically observed manifestation of hypertrophic cardiomyopathy in patients rather than TAC-induced hypertrophy. Furthermore, there are FIG. 7. Kurtosis along the principal eigenvectors (v 1 , v 2 , and v 3 ) of the diffusion tensor, averaged over the ROIs for healthy, sham, and TAC hearts. Kurtosis was computed from the beta distribution model. Table 2 Selected Parameter Values in Healthy, Sham, and TAC Hearts differences in species, tissue viability, contractility, fixation, and pulse sequence. Further investigation of such differences is needed for more direct comparison of data. Likewise, the community would benefit from greater harmonization of methods and reproducibility and validation studies.
Hypertrophic hearts exhibited significantly greater non-Gaussianity in the sheetlet and sheetlet-normal directions than healthy hearts. Specifically, an increase in low diffusivities was observed compared with the healthy and sham hearts. This may be caused by mesoscale tissue remodeling in response to pressure overload, including cardiomyocyte hypertrophy, decreased capillary density, and interstitial fibrosis (43).
The optimal model to describe the data was found to be highly dependent on maximum b-value. At b-values up to approximately 2000 s/mm 2 , the diffusion tensor model yielded the highest relative likelihood. As the maximum b-value increased, the effects of restricted diffusion are increasingly affecting measured data, and the truncated Oðb 2 Þ terms begin to become nonnegligible. Thus, the non-Gaussian models start to fit the data better. At the highest b-value considered in this study (10,000 s/mm 2 ), the beta distribution model offered the best fit, with slightly lower RMSE values than the biexponential model. Although we do not suggest that the distribution of diffusivities necessarily matches either of the distributions suggested by the models-it is clear from Figure 1 that vastly different distributions can yield almost identical signal attenuation profiles-we note that the general trends of these models are in good agreement. Neither model suggests that the maximum diffusivity in the TAC hearts was lower than in healthy or sham hearts, but rather that there was a greater volume fraction of low diffusivities perpendicular to locally prevailing cardiomyocyte orientations. The biexponential model was frequently the best model in regions exhibiting partial volume effects, most likely because this model has the flexibility to represent two distinct populations of diffusion. Although the diffusion kurtosis model frequently did not yield the highest relative likelihood, it was useful for providing initial values for other models, because it can be fit using linear least squares regression. Overall, the truncated Gaussian model offered the best fit out of the 13-parameter models.
In this study, models of non-Gaussian diffusion were fit in a diagonalized reference frame, based on the assumption that the principal directions of non-Gaussian diffusion are the same as those derived from the DTI model. Fitting in the diagonalized reference frame was necessary to make use of statistical models, which are only defined in a single dimension. Although our preliminary work (24) supports the assumption that Gaussian and non-Gaussian diffusion are well aligned in the myocardium, this removes the possibility to derive the kurtosis tensor (8) and therefore to estimate kurtosis in arbitrary orientations. It also may limit our understanding of the non-Gaussian behavior of regions with multiple populations of differently oriented cells, such as at the right ventricular wall or papillary muscle insertion points.
In this study, gadolinium chelate was used to reduce tissue T 1 and scan times. Gadolinium chelate is a known extracellular contrast agent in vivo, and differential localization in the intra-and extracellular compartments may affect the diffusion MRI signal. However, the distribution of gadolinium chelate ex vivo depends on the integrity of cell membranes after fixation. One in vitro study found that the T 1 of fixed cells doped with gadolinium chelate was reduced by >20 Â relative to fresh cells doped with gadolinium chelate and suggested that this finding constituted evidence of structural and functional alterations to the cell and nuclear membranes with formaldehyde fixation (44). Another study reported that the transmembrane water exchange rate was dramatically increased after fixation with Karnovsky's fixative (45), as we have used. These suggest that gadolinium chelate may have diffused into the intracellular space, thereby shortening T 1 in both intra-and extracellular compartments. It is also possible that with sufficiently high water exchange rates and long diffusion times, the intracellular and extracellular water may become well mixed. Further investigation would be needed for verification.
At present, it is technically challenging to measure non-Gaussian diffusion in the heart in vivo. In clinical diffusion imaging, non-Gaussian models have the potential to help in assessing myocardial perfusion (19). However, gradient hardware in clinical scanners typically limit b-values to 500 s/mm 2 (46)(47)(48), reducing sensitivity to non-Gaussian diffusion arising from finer tissue structures as demonstrated here. Stimulated echo approaches afford higher b-values for a given maximum gradient amplitude due to long diffusion times. Longer D would increase the interactions of water molecules with cellular restrictions and likely enhance non-Gaussian effects, particularly in the directions of the second and third eigenvectors. However, stimulated echo sequences are subject to myocardial strain, necessitating strain correction (48) or limiting imaging to temporal "sweet spots" (49). Higher b-values result in lower signal, thus higher SNR is needed to avoid noise floor bias. Furthermore, one DW readout is typically acquired at every other cardiac cycle, and a large number of b-values are required for non-Gaussian parameter estimation, increasing acquisition time. These drive the need for better SNR efficiency. Ongoing developments in high performance gradient systems, and the use of spin echo methods (50), convex optimized diffusion encoding (51), simultaneous multislice imaging (52), and compressed sensing (34) promise to improve SNR efficiency and therefore feasibility of in vivo cardiac non-Gaussian diffusion imaging.

CONCLUSION
This study represents the first systematic study of non-Gaussian diffusion in cardiac tissue. Several models of diffusion were presented and applied to fixed ex vivo rat hearts. At b-values up to approximately 2000 s/mm 2 , the diffusion tensor model yielded the highest relative likelihood, but reflects the signal poorly in the high b-value regime. The best model at high b-values, as selected by the AIC, was the beta distribution model. Differences in kurtosis and skewness, but not in DTI-derived parameters, were observed between healthy and hypertrophic hearts, demonstrating the higher sensitivity of the beta distribution model to pathology compared with the diffusion tensor model. We conclude that non-Gaussian models best represent the diffusion MRI signal in the myocardium at b-values exceeding 2000 s/mm 2 , and that parameters derived from these models may serve as useful biomarkers for assessing cardiac disease and remodeling.

APPENDIX
Here we provide additional details about a number of mathematical functions described in the Theory section.
The error function is related to the integral of the Gaussian distribution, and is defined as follows: The gamma function is an extension of the factorial function from positive integers to real and complex numbers. It is also known as the Euler integral of the second kind and is given as follows: x tÀ1 e Àx dx: The beta function is the Euler integral of the first kind and is given by Bðx; yÞ ¼ Z 1 0 t xÀ1 ð1 À tÞ yÀ1 dt ¼ GðxÞGðyÞ Gðx þ yÞ : [A3] Kummer's function is a standard form of confluent hypergeometric function and is given by M ða; b; zÞ ¼ X 1 n¼0 a ðnÞ z n b ðnÞ n! ; [A4] where a ð0Þ ¼ 1, a ðnÞ ¼ aða þ 1Þða þ 2Þ Á ÁÁða þ n À 1Þ denotes the rising factorial.