 Technical notes
 Open Access
 Published:
Bayesian intravoxel incoherent motion parameter mapping in the human heart
Journal of Cardiovascular Magnetic Resonance volume 19, Article number: 85 (2017)
Abstract
Background
Intravoxel incoherent motion (IVIM) imaging of diffusion and perfusion in the heart suffers from high parameter estimation error. The purpose of this work is to improve cardiac IVIM parameter mapping using Bayesian inference.
Methods
A secondorder motioncompensated diffusion weighted spinecho sequence with navigatorbased slice tracking was implemented to collect cardiac IVIM data in early systole in eight healthy subjects on a clinical 1.5 T CMR system. IVIM data were encoded along six gradient optimized directions with bvalues of 0–300 s/mm^{2}. Subjects were scanned twice in two scan sessions one week apart to assess intrasubject reproducibility. Bayesian shrinkage prior (BSP) inference was implemented to determine IVIM parameters (diffusion D, perfusion fraction F and pseudodiffusion D*). Results were compared to leastsquares (LSQ) parameter estimation. Signaltonoise ratio (SNR) requirements for a given fitting error were assessed for the two methods using simulated data. Reproducibility analysis of parameter estimation invivo using BSP and LSQ was performed.
Results
BSP resulted in reduced SNR requirements when compared to LSQ in simulations. Invivo, BSP analysis yielded IVIM parameter maps with smaller intramyocardial variability and higher estimation certainty relative to LSQ. Mean IVIM parameter estimates in eight healthy subjects were (LSQ/BSP): 1.63 ± 0.28/1.51 ± 0.14·10^{−3} mm^{2}/s for D, 13.13 ± 19.81/13.11 ± 5.95% for F and 201.45 ± 313.23/13.11 ± 14.53·10^{−3} mm^{2}/s for D ^{∗}. Parameter variation across all volunteers and measurements was lower with BSP compared to LSQ (coefficient of variation BSP vs. LSQ: 9% vs. 17% for D, 45% vs. 151% for F and 111% vs. 155% for D ^{∗}). In addition, reproducibility of the IVIM parameter estimates was higher with BSP compared to LSQ (BlandAltman coefficients of repeatability BSP vs. LSQ: 0.21 vs. 0.26·10^{−3} mm^{2}/s for D, 5.55 vs. 6.91% for F and 15.06 vs. 422.80·10^{−3} mm^{2}/s for D*).
Conclusion
Robust freebreathing cardiac IVIM data acquisition in early systole is possible with the proposed method. BSP analysis yields improved IVIM parameter maps relative to conventional LSQ fitting with fewer outliers, improved estimation certainty and higher reproducibility. IVIM parameter mapping holds promise for myocardial perfusion measurements without the need for contrast agents.
Background
Cardiovascular magnetic resonance (CMR) diffusion weighted imaging relies on signal attenuation due to random motion of water molecules in the presence of diffusion encoding gradients. Additionally, microvascular perfusion can contribute to the signal loss as described by the intravoxel incoherent motion (IVIM) model [1, 2]. According to Le Bihan et al. [3, 4], perfusion can be modeled as pseudo diffusion on a macroscopic scale, assuming random orientation of microvasculature in the capillary network. Consequently, the signal intensity can be described by a biexponential signal decay as a function of the diffusion encoding strength (bvalue). As the IVIM method is an endogenous contrast technique, its application is particularly suited to obtain a tissue perfusion surrogate where contrast agent administration is contraindicated. In recent years, this technique has gained significant momentum with successful applications in various body parts [5,6,7,8,9]. Beyond body and brain applications, IVIM of the invivo human heart has also been demonstrated [10, 11].
Cardiac IVIM may allow to delineate infarcted and ischemic areas showing good agreement with lategadolinium enhanced imaging [12, 13]. Moreover, IVIM may enable the assessment of chronic and acute ischemia [14] as well as conditions related to microvascular obstruction of the myocardium [15].
Despite recent progress, invivo cardiac diffusion weighted imaging still remains challenging due to cardiac and respiratory motion. Additionally, low signaltonoise ratio (SNR) and long scan times are major impediments to a wider acceptance in a clinical setting. Motion induced signal loss in spinecho (SE) based cardiac diffusion weighted imaging has been addressed by firstorder motioncompensated diffusion gradient designs in conjunction with careful cardiac trigger delay selection [16] and more recently by secondorder motion compensation [17,18,19]. Initial results of the application of secondorder motion compensation for IVIM acquisitions during systole have previously been presented in a porcine model [14]. For diastolic imaging, timeshifted triggering and dedicated post processing using principal component analysis (PCA) filtering in combination with temporal maximum intensity projection (PCATMIP) has been proposed [20,21,22,23].
Experimentally, cardiac IVIM parameters were initially reported for the invivo canine heart by Callot et al. [24]. The measured diffusion weighted signal agreed well with the biexponential IVIM model with reduced signal decay in the absence of perfusion postmortem [14].
IVIM parameter maps of various organs such as brain or heart [20, 25] are typically of noisy appearance. Due to the nonlinearity and bad conditioning of the regression problem, the perfusion related parameters are estimated with considerable error at typical SNR values as shown in [26, 27]. Besides modifying the data acquisition protocol to obtain higher SNR at the expense of lower spatial resolution and/or longer scan time, group analysis of longitudinal data of individuals incorporating both intra and intersubject variations [28] or regional smoothing [29] have been proposed. These approaches are, however, limited by the necessity of repeated measurements across multiple independent subjects or loss of spatial resolution and increased partialvoluming, respectively.
To address the SNR limitation of IVIM analysis, a hierarchical Bayesian data analysis framework has been presented by Orton et al. [30] and demonstrated for liver application. Using this approach, information across the regionofinterest is taken into account for voxelwise parameter inference. Parameter estimation is performed using a posterior distribution combining data likelihood and a hierarchical prior. This combination enables effective denoising of parameter maps with reduced parameter estimation error.
The objective of the present work was to implement and assess Bayesian shrinkage prior (BSP) inference for IVIM parameter mapping of the invivo human heart and compare its performance to segmented leastsquares (LSQ) fitting.
Theory
Intravoxel incoherent motion
The IVIM model [1] in Eq. (1) describes signal magnitude of diffusion weighted images as biexponential decay. In addition to diffusion induced signal attenuation, a second compartment of perfusion induced pseudodiffusion is taken into account:
where S(b) describes the measured signal as a function of bvalue, S _{0} the signal without diffusion weighting (b=0 s/mm^{2}), D the diffusion constant, F the perfusion fraction and D ^{∗} the pseudodiffusion constant. Note that capital F is used for the perfusion fraction to be consistent with the notation of Orton et al. [30].
Leastsquares fitting
For LSQ fitting, a segmented approach [6] is implemented assuming the contribution of the perfusion to reach a maximum of F/(1 − F) at b=0 s/mm^{2} and to drop to negligible values for bvalues b ≫ b _{ Split }. In practice, high bvalues (b ≥ b _{ Split }=200 s/mm^{2}) are fitted to a monoexponential diffusiononly model:
If a nondiffusion weighted image S _{0} is not available, the intercept S _{ int } does not allow for a direct calculation of the perfusion fraction F as described in [6], but S _{ int } = S _{0} ∙ (1 − F) enables to eliminate S _{0} in the biexponential model. In the second step of the segmented regression, the perfusion related parameters F and D ^{∗} are estimated using the predetermined diffusion coefficient D and the monoexponential intercept S _{ int } while taking into account all considered bvalues. By substituting S _{0} in Eq. (1), the signal model reads accordingly:
The nonlinear regression is implemented using an interiorpoint algorithm in Matlab (Mathworks, Natick, Massachusetts, USA) and constrained by an inequality together with box constraints as in [30]:
Bayesian shrinkage prior inference
For Bayesian inference as presented in [30], a marginalized data likelihood is used along with a multivariate Gaussian prior combined with Jeffrey’s prior [31]. The approach is implemented using a Markov chain Monte Carlo (MCMC) method as described in the Appendix 1.
Methods
Computer simulations
IVIM parameter ranges were simulated according to values reported for the invivo heart [14, 20]. The diffusion coefficient was set to D=1.5·10^{−3} mm^{2}/s, while three perfusion regimes (low, intermediate, high) were considered (F/D ^{∗}=10/10, 15/15, 20/20 %/10^{−3} mm^{2}/s). The simulated SNRs ranged from 10 to 100 in steps of 10 and from 100 to 200 in steps of 25. Gaussian distributed noise was added followed by magnitude detection to yield Rician distributed noise mimicking the noise distribution of CMR magnitude images. A single Monte Carlo simulation run consisted of 1,000 IVIM data sets with bvalues as used in the invivo part of this study: 20 to 100 s/mm^{2} in steps of 20 s/mm^{2}, 125 to 200 s/mm^{2} in steps of 25 s/mm^{2}, 250 and 300 s/mm^{2}. Both bias \( \left\left\langle \widehat{p}\right\rangle {p}_{Ref}\right/{p}_{Ref} \) and variation \( {\widehat{\sigma}}_p/{p}_{Ref} \) with p = D, F and D ^{∗}, p _{Ref} the simulated parameter, \( \left\langle \widehat{p}\right\rangle \) the mean estimate and \( {\widehat{\sigma}}_p \) the standard deviation of the estimated parameters were calculated and are reported as relative errors. The simulation was repeated and resulting parameter estimation errors were averaged 100 times.
Invivo measurements
Secondorder motioncompensated SE diffusion weighted imaging [17, 18] was implemented on a 1.5 T CMR system (Achieva, Philips Healthcare, Best, The Netherlands), see Fig. 1. Signal was received with a 5 channel cardiac receiver array. Written informed consent was obtained from all subjects prior to imaging. The study protocol was approved by the ethics committee of the Canton of Zurich. Consent included imaging as well as publication of anonymized data.
Data were acquired in eight healthy subjects without history of cardiac disease (6 female, 2 male, weight 64 ± 8 kg, age 26 ± 4 years, heart rate 64 ± 9 beats/min, min/max heart rates: 51/81 beats/min) on two separate occasions. Prior to diffusion imaging, cine data with a temporal resolution of 10 ms were acquired in two chamber and short axis view orientations. Using the cine images, systolic quiescent time points were determined visually on a per subject basis.
CMR diffusion weighted imaging was performed during freebreathing in shortaxis view orientation using singleshot EPI readout with the reduced fieldofview (FOV) technique locallook (LoLo) [32]. Slice tracking to account for breathing motion was controlled by a respiratory 1D navigator pencil beam placed on the right hemi diaphragm, accepting all data. A 1–3–31 binomial spectralspatial excitation pulse for fat suppression [33] was employed. Images were acquired with inplane resolution: 2.4 × 2.4 mm^{2}, slice thickness: 10 mm, one midventricular slice, FOV: 230 × 105 mm^{2}, acquired kspace lines: 43, TR/TE: 2 RR/73 ms, flip angle: 81 ± 1° (heart rate dependent Ernst angle [34]), 8 signal averages and 6 vendor gradient optimized diffusion encoding directions. The applied diffusion encoding strengths included the values described in the previous sections (20 – 100 s/mm^{2} in steps of 20 s/mm^{2}, 125 – 200 s/mm^{2} in steps of 25 s/mm^{2}, 250 and 300 s/mm^{2}) together with 0 s/mm^{2}. The trigger delay for the SE sequence was set to 25% peak systolic contraction [17] with a mean trigger delay of 78 ± 3 ms. Acquisition of the 8 signal averages for each diffusion encoding strength and direction was equally distributed along the measurement. Total scan time was about 18 min at a heart rate of 60 beats per minute.
Invivo SNR measurements were performed in each subject. To measure noise, the scans were repeated with radiofrequency and gradient pulses switched off. Sufficient time (>10 s) was allowed between image and noise acquisition to ensure complete signal decay. SNR was determined for each voxel as described in [35].
Imaging in each subject was repeated in consecutive sessions separated by one week to assess intrasubject reproducibility.
In addition, diffusion data in an animal model of myocardial infarction was evaluated; see Appendix 2 for further details.
Data postprocessing
For invivo IVIM parameter mapping, images were first registered using a dedicated groupwise image registration method [36] employing total variation displacement regularization and a PCAbased image similarity metric [37] to correct for inplane residual geometric inconsistencies. Afterwards, complex averaging [38] of the signal averages was performed. The IVIM parameters of both regression methods (LSQ and BSP) and SNR were determined upon manual masking the left ventricular myocardium. The same segmentation was used for both regression methods (LSQ and BSP). To avoid partial volume effects, voxels at the epi and endocardial borders were excluded during the segmentation process and all voxels within the segmented regionofinterest were used for further analysis. The image magnitudes were corrected for heart rate variations using recorded RR intervals and published T_{1} values of the myocardium [39]. Figure 2 summarizes all post processing steps. IVIM analysis was performed on data with b ≥ 20 s/mm^{2} to suppress artifacts from blood flow while mean SNR was determined on b=0 s/mm^{2} images.
For BSP inference, the total number of MCMC samples was set to N _{ s }=20,000. A “burnin” period of 10,000 (discarded) samples was used before actual sampling. The Markov chains were initialized with LSQ estimates of the IVIM parameters. Note that the Markov chains can be started from arbitrary starting values, however a starting point close to the actual parameter estimates shortens the burnin phase and hence saves computation time. Further details of the estimation method can be found in the Appendix of Orton et al. [30]. A vectorized approach of the referenced procedure was implemented in Matlab (Mathworks) and run on standard PC hardware (2.9 GHz, 16 GB RAM).
Reproducibility analysis
In order to assess reproducibility of two consecutive scans, BlandAltman analysis was performed and the coefficient of variability was calculated for both scan sessions.
Results
Computer simulations
In Fig. 3, relative errors of D, F and D ^{∗} for BSP versus LSQ as a function of SNR are reported for a Monte Carlo simulation. Both methods show overall decreasing errors for increasing SNR.
The bias of D with both LSQ and BSP is reduced to below 20% for all perfusion scenarios at SNR ≥ 20, with the bias of LSQ remaining between 3 and 10% even at high SNRs. The variation of D with BSP is consistently lower over the entire SNR range compared to LSQ. It drops below 20% error at a SNR of approx. 90 for LSQ, but remains below 20% already for the lowest simulated SNR of 10 for BSP. Estimation of the perfusion fraction F yields lower bias with BSP relative to LSQ between an SNR of 20 to 90–150 depending on the perfusion regime tested. LSQ shows an increase in bias for SNR ≥ 125–175. The variation of F with BSP is consistently lower compared to LSQ over the entire SNR range and perfusion regimes. The relative error is below 20% for SNR ≥ 30 if BSP is used for inference. Depending on the perfusion regime simulated, the error using LSQ remains above that threshold except for the high perfusion regime at a SNR ≥ 175. The SNR dependency of D ^{∗} shows consistently lower bias and variation for BSP relative to LSQ for SNR ≥ 20. For LSQ, bias remains above 20% error for an SNR of 200 for all perfusion regimes. Variation of D ^{∗} with LSQ remains also above 20% even at an SNR of 200 for all perfusion regimes while bias and variation with BSP are bound to below 20% for SNR values above 40 and 60, respectively. Based on the simulation, an overall minimum SNR of 30–60 depending on the perfusion regime is identified for BSP to determine D, F, and D ^{∗} within 20% bias and variation. The LSQ method exhibits errors above the mentioned threshold even at a SNR of 200.
Invivo measurements
The invivo SNR measured without diffusion weighting was 19 ± 3 for one signal average, resulting in an SNR of approximately 54 for averaged data. Figure 4 shows example invivo magnitude images of all bvalues. The bright blood pool signal in the center of the image is dephased with increasing diffusion weighting. Example trace magnitude signals averaged across the regionofinterest are displayed. In addition, trace signals from all volunteers and repetitions are plotted.
In Fig. 5, example IVIM parameter maps computed with LSQ and BSP along with corresponding histograms are shown. While LSQ maps exhibit spatial noise and patchlike structures, BSP yields a more uniform distribution in the myocardium which is reflected in narrower distributions of D, F and D ^{∗}. Of note, LSQ resulted in a high number of voxels in which the estimated IVIM parameters reached or were close to the box constraints.
Figure 6 summarizes various parameter estimates together with regression quality measures (LSQ red boxes, BSP blue boxes) as in [30]. The left column summarizes mean and median estimates across the corresponding regionsofinterest of all parameters. Both the LSQ mean and median estimates of D tend to relatively high values compared to BSP, while the prior mean of BSP is within the range of previously reported values [20, 40]. Considering the parameter F, there are notable differences among the mean and median estimates, indicating the presence of fitting outliers. Again, the prior mean of BSP is in the range of previously reported values [20]. The mean estimates of D ^{∗} are strongly influenced by (high valued) outliers in the regionofinterest, explaining the difference between mean and median LSQ estimates. The BSP prior means take values close to the LSQ median values. The variability measures within the regionofinterest in the middle column show reduced variability for BSP versus LSQ in all parameters both considering standard deviation and percentile based measures. The fit quality in terms of median estimated standard deviation under the posterior (for example σ _{ d }) is displayed in the right column. The BSP based deviations are consistently lower compared to the LSQ based values.
For reproducibility analysis, medians across the left ventricular myocardium/regionofinterest were considered because of the large amount of outliers for LSQ fitting. Figure 7 shows the BlandAltman analysis of two consecutive scans within one session. Mean biases (LSQ/BSP) of −0.02/−0.05·10^{−3} mm^{2}/s for D, −0.58/+0.51% for F and +26.28/−1.56·10^{−3} mm^{2}/s for D ^{∗} were found. The BlandAltman coefficients of repeatability are (LSQ/BSP): 0.26/0.21·10^{−3} mm^{2}/s for D, 6.91/5.55% for F and 422.80/15.06·10^{−3} mm^{2}/s for D*.
Figure 8 shows a summary over all measurement estimates. The upper row displays medians across the regionsofinterest for both sessions. As in Fig. 6, the estimates of D were found to be higher for LSQ compared to BSP. The medians of the LSQ/BSP estimates are covering ranges of 0.61/0.51·10^{−3} mm^{2}/s for D, 14.79/10.27% for F, 763.37/27.42·10^{−3} mm^{2}/s for D ^{∗}. The lower row of Fig. 8 reports all measurements by displaying the means across all voxels within the regionsofinterest and the corresponding standard deviations. The mean values of the IVIM parameters are (LSQ/BSP) 1.63 ± 0.28/1.51 ± 0.14·10^{−3} mm^{2}/s for D, 13.13 ± 19.81/13.11 ± 5.95% for F and 201.45 ± 313.23/13.11 ± 14.53·10^{−3} mm^{2}/s for D ^{∗}. The estimates for all IVIM parameters from the two inference procedures are significantly different (p < 0.05) from each other using the Wilcoxon signedrank test. Both mean D and F are within 10% relative difference, but the mean estimates of D ^{∗} are one order of magnitude different from each other. This is again due to the high number of outliers produced by the LSQ method. The standard deviations of all three IVIM parameters are consistently lower for BSP compared to LSQ.
Potential scan time reduction was investigated by retrospectively skipping diffusion encoding gradient directions. The reduced SNR due to data subsampling leads to a mean absolute error across all voxels of all measurements of 0.25/0.14·10^{−3} mm^{2}/s (D), 11.38/4.67% (F) and 185.38/15.82·10^{−3} mm^{2}/s (D ^{∗}) for LSQ/BSP if only three instead of all six directions are used as shown in Fig. 9.
The LSQ and BSP estimates were further compared by a regression analysis of the median estimates across the corresponding regionsofinterest as shown in Fig. 10. The diffusion coefficient exhibits an approximately linear correlation between the two methods with LSQ estimates tending to higher values. For the perfusion parameters there is no clear linear correlation. Especially D ^{∗} exhibits outliers. The KullbackLeibler divergences D _{ KL }(LSQ  BSP) of the BSP parameter estimates from the LSQ estimates are summarized for all measurements in Fig. 10. The median divergences and standard deviations were: 15 ± 11/16 ± 11/14 ± 10 bit for D/F/D ^{∗} respectively.
The infarcted septal region in an animal model exhibits reduced blood flow in a conventional contrast enhanced first pass perfusion scan as well as reduced IVIM perfusion parameters for both regression methods. The BSP derived maps do not contain outliers and hence allow a clearer delineation of the infarcted area compared to LSQ. Further details can be found in Appendix 2.
Discussion
In the present work, Bayesian shrinkage prior inference has been implemented and compared to segmented leastsquares fitting for IVIM parameter mapping in the invivo human heart.
Robust data acquisition was possible using a secondorder motioncompensated diffusion weighted spinecho sequence [17] triggered to early systole. Using a trigger delay of 25% peak systole is advantageous because of increased coronary flow compared to peak contraction. Moreover, a large part of systole is potentially available for IVIM acquisitions as shown in [17]: trigger delays in the range of 15–77/79% peak systole at the apex/base allow for robust diffusion data acquisition. In addition, imaging in systole has the advantage of a relatively thick myocardium compared to the voxel size.
Using motioncompensated diffusion gradients may lead to a reduced sensitivity to blood circulation and myocardial perfusion [41]. However, deviations from a monoexponential diffusion model at lower bvalues due to perfusion were observed in all measurements indicating sufficient sensitivity to microcirculation and perfusion. Balancing motioninduced signal loss due to cardiac bulk motion while achieving optimal sensitivity to perfusion is however a subject deserving further attention.
Computer simulations revealed minimum SNR thresholds of 30–60 for relative errors in terms of bias and variation of 20% each, depending on the perfusion regime for BSP while the LSQ method required a minimum SNR of at least 200. The increase in bias in F for the SNR range of 125–175 using LSQ is suspected to be an artefact of the segmented fit, which leads to an error propagation of D and S _{ int } into the estimation of the perfusion parameters. Furthermore, we note that Federau and colleagues also reported a similar increase in bias in F in Fig. 1 of [27] using a segmented approach; albeit in the SNR range of approx. 20–40 with simulated IVIM parameters which are commonly found in the brain (D=0.7·10^{−3}mm^{2}s, F=4% and D ^{∗}=17·10^{−3}mm^{2}s). Even though the relative bias in F in the considered SNR range is below about 20%, the segmented fit might benefit from a joint parameter estimation (potentially using two regimes of mono and biexponential decay) in this regard. While these simulation results are indicative, several factors confounding invivo measurements have not been taken into account in the simulations including residual motion artifacts and partial voluming with hyperintense blood signal and epicardial fat [42]. These effects would lead to a broadening of the parameter histograms. Accordingly, the width of the prior is increased by the presence of a large number of affected voxels. The shrinking procedure in these cases is less effective.
Invivo, BSP analysis resulted in IVIM parameter maps with considerably smaller intramyocardial standard deviations relative to LSQ. Both variability and estimation uncertainty in terms of standard deviation under the posterior were greatly reduced with BSP compared to LSQ (Fig. 6), indicating the benefit of taking into account prior knowledge. The setting of arbitrary fit constraints was obsolete in the BSP procedure. In addition, BSP regression was aided by the prior which led to the elimination of outliers on or close to the fit boundaries.
An effective spatial denoising of the parameters can be achieved because the prior in BSP is chosen to be a unimodal distribution. This prior assumes a population mean of the IVIM parameters for the whole regionofinterest and hence assumes the myocardium of the left ventricle to have a rather homogenous spatial tissue characteristic. If local pathologies such as myocardial infarcts (as for example presented in Appendix 2) and corresponding fibrous tissue with reduced perfusion [14] are present, a different choice of priors such as a spatial homogeneous prior [43] is deemed more appropriate [42]. Alternatively, multimodal priors can be applied to distinguish among different tissue types while still retaining spatial information. Methods using mixture models [44] of multivariate Gaussians could also be implemented to address this limitation. In contrast to the LSQ approach, which allows data processing on a pixelbypixel basis, the BSP method requires presegmentation of the data, which renders automation in a postprocessing workflow more challenging.
Overall, the invivo IVIM parameters measured in this study are in good accordance with recent literature [10, 20]. The diffusion coefficients D found in this study using LSQ and BSP were within the range of the values found in [10, 20]. A higher measured diffusivity is indicative of residual motion effects in the data [45]. In addition, the mean diffusivity measured using spinecho based diffusion tensor imaging (DTI) during early systole [40] was within 14% and 6% of the mean measured diffusion coefficient of the present study for LSQ and BSP, respectively. The measured perfusion fraction F was found to be lower compared to data reported in [10]. The LSQ and BSP estimation of approx. 13% is close to the 12% found in [20] during diastole. The pseudodiffusion coefficients D ^{∗} of 201.50 ± 313.20 mm^{2}/s (LSQ) and 13.11 ± 14.53 mm^{2}/s (BSP) as measured in the present study are different from previous data (52.68 ± 52.61·10^{−3} mm^{2}/s [10], 43.6 ± 9.2·10^{−3} mm^{2}/s [20]). However, D ^{∗} usually contains the highest number of outliers and the mean across the regionofinterest is therefore prone to be heavily influenced by the choice of the actual value of the boxconstraints. If the medians across the regionsofinterest are considered, it is shown that the majority of the LSQ estimates gather in the range of 20 to 50·10^{−3} mm^{2}/s. Moreover, data in Fig. 6 indicates that mean estimates of D ^{∗} are considerably influenced by outliers with median parameter values close to the ones found using the BSP method.
The IVIM perfusion parameters in an infarcted animal model (see also Appendix 2) show good accordance with the perfusion defect visible in the contrast enhanced perfusion scan, especially considering the blood flow related [25, 27] product of F × D ^{∗} of the BSP derived estimates. However, the extension of the infarct indicated by IVIM appears smaller compared to the darkened area of the contrast enhanced perfusion scan. This might be due to residual motion and/or partial voluming which can yield elevated IVIM parameter estimates.
All SNR measurements were obtained from a single signal average. Thereby confounding factors due to image registration and phase correction for averaging of complex data were avoided. The invivo SNR was above the 20% parameter error threshold found in simulations. The scan time of ca. 18 min for a heart rate of 60 bpm for this study was inbetween previous scanning times of 15 min [10] and 20 min [20]. Optimizations of the experiment design in terms of bvalue distribution [46], higher static field strength and improved gradient performance to reduce echo times may allow to reduce parameter estimation error and scan time.
By design, Bayesian approaches exploiting information across the regionofinterest can be used to examine distributed rather than focal pathologies of the myocardium. Accordingly, potential applications relate to microvascular obstruction and reduced and/or delayed perfusion of the myocardium in hypertrophic cardiomyopathy and diabetes [15, 47].
Conclusion
Bayesian IVIM parameter mapping yields improved parameter maps relative to conventional segmented leastsquares fitting in the human heart. In conjunction with motioncompensated diffusion weighted spinecho sequences, robust parameter estimation can be achieved providing a tissue perfusion surrogate without contrast agent application. Further invivo studies are now warranted to assess the performance of the method in relevant patient populations.
Abbreviations
 BSP:

Bayesian shrinkage prior
 FOV:

Fieldofview
 IVIM:

Intravoxel incoherent motion
 LoLo:

Locallook
 LSQ:

Leastsquares
 MCMC:

Markov chain Monte Carlo
 PCA:

Principal component analysis
 SNR:

Signaltonoise ratio
 TE:

Echo time
 TR:

Repetition time
References
 1.
Le Bihan D. Intravoxel incoherent motion imaging using steadystate free precession. Magn Reson Med. 1988;7:346–51.
 2.
Le Bihan D. Intravoxel incoherent motion perfusion MR imaging: a wakeup call. Radiology. 2008;249:748–52. https://doi.org/10.1148/radiol.2493081301.
 3.
Le Bihan D, Breton E, Lallemand D, Aubin ML, Vignaud J, LavalJeantet M. Separation of diffusion and perfusion in intravoxel incoherent motion MR imaging. Radiology. 1988;168:497–505. https://doi.org/10.1148/radiology.168.2.3393671.
 4.
Le Bihan D, Turner R. The capillary network: a link between IVIM and classical perfusion. Magn Reson Med. 1992;27:171–8.
 5.
Yamada I, Aung W, Himeno Y, Nakagawa T, Shibuya H. Diffusion coefficients in abdominal organs and hepatic lesions: evaluation with intravoxel incoherent motion echoplanar MR imaging. Radiology. 1999;210:617–23. https://doi.org/10.1148/radiology.210.3.r99fe17617.
 6.
Notohamiprodjo M, Chandarana H, Mikheev A, Rusinek H, Grinstead J, Feiweier T, Raya JG, Lee VS, Sigmund EE. Combined intravoxel incoherent motion and diffusion tensor imaging of renal diffusion and flow anisotropy. Magn Reson Med. 2015;73:1526–32. https://doi.org/10.1002/mrm.25245.
 7.
Lemke A, Laun FB, Klau M, Re TJ, Simon D, Delorme S, Schad LR, Stieltjes B. Differentiation of pancreas carcinoma from healthy pancreatic tissue using multiple bvalues. Investig Radiol. 2009;44:769–75. https://doi.org/10.1097/RLI.0b013e3181b62271.
 8.
Luciani A, Vignaud A, Cavet M, et al. Liver cirrhosis: intravoxel incoherent motion MR imaging—pilot study. Radiology. 2008;249:891–9. https://doi.org/10.1148/radiol.2493080080.
 9.
Patel J, Sigmund EE, Rusinek H, Oei M, Babb JS, Taouli B. Diagnosis of cirrhosis with intravoxel incoherent motion diffusion MRI and dynamic contrastenhanced MRI alone and in combination: preliminary experience. J Magn Reson Imaging. 2010;31:589–600. https://doi.org/10.1002/jmri.22081.
 10.
Froeling M, Strijkers GJ, Nederveen AJ, Chamuleau SA, Luijten PR. Feasibility of in vivo whole heart DTI and IVIM with a 15 minute acquisition protocol. J Cardiovasc Magn Reson. 2014;16:O15. https://doi.org/10.1186/1532429X16S1O15.
 11.
Moulin KK, Croisille P, Feiweier T, Delattre BMA, Wei H, Robert B, Beuf O, Viallon M. In vivo freebreathing DTI and IVIM of the whole human heart using a realtime slicefollowed SEEPI navigatorbased sequence: a reproducibility study in healthy volunteers. Magn Reson Med. 2016;76:70–82. https://doi.org/10.1002/mrm.25852.
 12.
Deux JF, Maatouk M, Vignaud A, Luciani A, Lenczner G, Mayer J, Lim P, DuboisRandé JL, Kobeiter H, Rahmouni A. Diffusionweighted echo planar imaging in patients with recent myocardial infarction. Eur Radiol. 2011;21:46–53. https://doi.org/10.1007/s0033001019126.
 13.
Laissy JP, Gaxotte V, Irondelaissy E, Klein I, Ribet A, Bendriss A, Chillon S, SchoumanClaeys E, Steg PG, Serfaty JM. Cardiac diffusionweighted MR imaging in recent, subacute, and chronic myocardial infarction: a pilot study. J Magn Reson Imaging. 2013;38:1377–87. https://doi.org/10.1002/jmri.24125.
 14.
Deuster C Von, Stoeck CT, Wissmann L, Spinner G, Fleischmann T, Emmert MY, Cesarovic N, Kozerke S. Verification of the intravoxel incoherent motion (IVIM) model in the porcine heart. In: Proc Intl Soc Mag Reson Med. 2015;23.
 15.
Ismail TF, Hsu LY, Greve AM, et al. Coronary microvascular ischemia in hypertrophic cardiomyopathy  a pixelwise quantitative cardiovascular magnetic resonance perfusion study. J Cardiovasc Magn Reson. 2014;16:1–10. https://doi.org/10.1186/s1296801400491.
 16.
Gamper U, Boesiger P, Kozerke S. Diffusion imaging of the in vivo heart using spin echoesconsiderations on bulk motion sensitivity. Magn Reson Med. 2007;57:331–7. https://doi.org/10.1002/mrm.21127.
 17.
Stoeck CT, von Deuster C, Genet M, Atkinson D, Kozerke S. Secondorder motioncompensated spin echo diffusion tensor imaging of the human heart. Magn Reson Med. 2016;75:1669–76. https://doi.org/10.1002/mrm.25784.
 18.
Welsh CL, DiBella EVR, Hsu EW. Higherorder motioncompensation for in vivo cardiac diffusion tensor imaging in rats. IEEE Trans Med Imaging. 2015;34:1843–53. https://doi.org/10.1109/TMI.2015.2411571.
 19.
Nguyen C, Fan Z, Sharif B, He Y, Dharmakumar R, Berman DS, Li D. In vivo threedimensional high resolution cardiac diffusionweighted MRI: a motion compensated diffusionprepared balanced steadystate free precession approach. Magn Reson Med. 2014;72:1257–67. https://doi.org/10.1002/mrm.25038.
 20.
Moulin KK, Croisille P, Feiweier T, Delattre BMA, Wei H, Robert B, Beuf O, Viallon M. In vivo freebreathing DTI and IVIM of the whole human heart using a realtime slicefollowed SEEPI navigatorbased sequence: a reproducibility study in healthy volunteers. Magn Reson Med. 2015;3:70–82. https://doi.org/10.1002/mrm.25852.
 21.
Pai VM, Rapacchi S, Kellman P, Croisille P, Wen HPCATMIP. Enhancing signal intensity in diffusionweighted magnetic resonance imaging. Magn Reson Med. 2011;65:1611–9. https://doi.org/10.1002/mrm.22748.
 22.
Rapacchi S, Wen H, Viallon M, Grenier D, Kellman P, Croisille P, Pai VM. Low bvalue diffusionweighted cardiac magnetic resonance imaging: initial results in humans using an optimal timewindow imaging approach. Investig Radiol. 2011;46:751–8. https://doi.org/10.1097/RLI.0b013e31822438e8.
 23.
Delattre BM, Viallon M, Wei H, Zhu Y, Pai VM, Wen H, Croisille P. Intravoxel incoherent motion applied to cardiac diffusion weighted MRI using breathhold acquisitions in healthy volunteers. J Cardiovasc Magn Reson. 2012;14:P261. https://doi.org/10.1186/1532429X14S1P261.
 24.
Callot V, Bennett E, Decking UKM, Balaban RS, Wen H. Vivo study of microcirculation in canine myocardium using the IVIM method. Magn Reson Med. 2003;50:531–40. https://doi.org/10.1002/mrm.10568.
 25.
Federau C, O’Brien K, Meuli R, Hagmann P, Maeder P. Measuring brain perfusion with intravoxel incoherent motion (IVIM): initial clinical experience. J Magn Reson Imaging. 2014;39:624–32. https://doi.org/10.1002/jmri.24195.
 26.
WC W, Chen YF, Tseng HM, Yang SC, My PC. Caveat of measuring perfusion indexes using intravoxel incoherent motion magnetic resonance imaging in the human brain. Eur Radiol. 2015;25:2485–92. https://doi.org/10.1007/s003300153655x.
 27.
Federau C, O’Brien K, Birbaumer A, et al. Functional mapping of the human visual cortex with intravoxel incoherent motion MRI Zochowski M. PLoS One. 2015;10:1–11. https://doi.org/10.1371/journal.pone.0117706.
 28.
Huang HM, Shih YY, Lin C. Formation of parametric images using mixedeffects models: a feasibility study. NMR Biomed. 2016;29:239–47. https://doi.org/10.1002/nbm.3453.
 29.
Rheinheimer S, Stieltjes B, Schneider F, Simon D, Pahernik S, Kauczor HU, Hallscheidt P. Investigation of renal lesions by diffusionweighted magnetic resonance imaging applying intravoxel incoherent motionderived parametersinitial experience. Eur J Radiol. 2012;81:e310–6. https://doi.org/10.1016/j.ejrad.2011.10.016.
 30.
Orton MR, Collins DJ, Koh DM, Leach MO. Improved intravoxel incoherent motion analysis of diffusion weighted imaging by data driven Bayesian modeling. Magn Reson Med. 2014;71:411–20. https://doi.org/10.1002/mrm.24649.
 31.
Jeffreys H. An invariant form for the prior probability in estimation problems. Proc R Soc Lond A Math Phys Sci. 1946;186:453–61.
 32.
Feinberg DA, Hoenninger JC, Crooks LE, Kaufman L, Watts JC, Arakawa M. Inner volume MR imaging: technical concepts and their application. Radiology. 1985;156:743–7. https://doi.org/10.1148/radiology.156.3.4023236.
 33.
Meyer CH, Pauly JM, Macovski A, Nishimura DG. Simultaneous spatial and spectral selective excitation. Magn Reson Med. 1990;15:287–304.
 34.
Ernst RR. Application of Fourier transform spectroscopy to magnetic resonance. Rev Sci Instrum. 1966;37:93–102. https://doi.org/10.1063/1.1719961.
 35.
NordmeyerMassner J A, De Zanche N, Pruessmann KP. Mechanically adjustable coil array for wrist MRI. Magn Reson Med. 2009;61:429–438. doi: https://doi.org/10.1002/mrm.21868.
 36.
Vishnevskiy V, Gass T, Szekely G, Tanner C, Goksel O. Isotropic total variation regularization of displacements in parametric image registration. IEEE Trans Med Imaging. 2016;62:1–1. https://doi.org/10.1109/TMI.2016.2610583.
 37.
Huizinga W, Poot DHJ, Guyader JM, et al. PCAbased groupwise image registration for quantitative MRI. Med Image Anal. 2015; https://doi.org/10.1016/j.media.2015.12.004.
 38.
Scott AD, NiellesVallespin S, Ferreira PF, McGill LA, Pennell DJ, Firmin DN. The effects of noise in cardiac diffusion tensor imaging and the benefits of averaging complex data. NMR Biomed. 2016;29:588–99. https://doi.org/10.1002/nbm.3500.
 39.
Messroghli DR, Plein S, Higgins DM, Walters K, Jones TR, Ridgway JP, Sivananthan MU. Human myocardium: singlebreathhold MR T1 mapping with high spatial resolutionreproducibility study. Radiology. 2006;238:1004–12. https://doi.org/10.1148/radiol.2382041903.
 40.
von Deuster C, Stoeck CT, Genet M, Atkinson D, Kozerke S. Spin echo versus stimulated echo diffusion tensor imaging of the in vivo human heart. Magn Reson Med. 2016;76:862–72. https://doi.org/10.1002/mrm.25998.
 41.
Wetscherek A, Stieltjes B, Laun FB. Flowcompensated intravoxel incoherent motion diffusion imaging. Magn Reson Med. 2015;74:410–9. https://doi.org/10.1002/mrm.25410.
 42.
While PT. A comparative simulation study of Bayesian fitting approaches to intravoxel incoherent motion modeling in diffusionweighted MRI. Magn Reson Med. 2017;77 https://doi.org/10.1002/mrm.26598.
 43.
Freiman M, PerezRossello JM, Callahan MJ, Voss SD, Ecklund K, Mulkern RV, Warfield SK. Reliable estimation of incoherent motion parametric maps from diffusionweighted MRI using fusion bootstrap moves. Med Image Anal. 2013;17:325–36. https://doi.org/10.1016/j.media.2012.12.001.
 44.
Stoneking CJ. Bayesian inference of Gaussian mixture models with noninformative priors; 2014. https://arxiv.org/abs/1405.4895.
 45.
Stoeck CT, von Deuster C, van Gorkum R, Kozerke S. Impact of eddycurrents and cardiac motion in DTI of the invivo heart  a comparison of secondorder motion compensated SE versus STEAM. In: Proc Intl Soc Mag Reson Med. 2017.
 46.
Lemke A, Stieltjes B, Schad LR, Laun FB. Toward an optimal distribution of bvalues for intravoxel incoherent motion imaging. Magn Reson Imaging. 2011;29:766–76. https://doi.org/10.1016/j.mri.2011.03.004.
 47.
Fang ZY, Prins JB, Marwick TH. Diabetic cardiomyopathy: evidence, mechanisms, and therapeutic implications. Endocr Rev. 2004;25:543–67. https://doi.org/10.1210/er.20030012.
 48.
Jones DK, Horsfield MA, Simmons A. Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging. Magn Reson Med. 1999;42:515–25.
Acknowledgements
The authors thank Dr. Sudhir Shankar Raman (Translational Neuromodeling Unit, University and ETH Zurich) for his help and advice on Bayesian statistics.
Funding
The authors acknowledge funding from the Seventh Framework Programme of the European Commission (VPHDARE@IT) and the Swiss National Science Foundation (320030_153014).
Availability of data and materials
The data that support the findings of this study are available from the corresponding author upon reasonable request subject to restriction on use by the Ethics Committee of the Canton of Zurich.
Author information
Affiliations
Contributions
SK, GS and CVD conceived of the study. CTS developed and implemented the diffusion weighting sequence. GS and CVD performed data acquisition. GS implemented the Bayesian inference procedure, performed simulations, processed data and drafted the manuscript. KT helped implementing the Bayesian inference procedure. All authors participated in revising the manuscript and read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
The study was approved by the Ethics Committee of the Canton of Zurich, Switzerland, and all subjects provided written informed consent. Imaging was performed at the Zurich University Hospital, Zurich, Switzerland. Anonymized data was analyzed at ETH Zurich with approval by the mentioned authority.
The animal experiment was performed in adherence to the Swiss law of animal protection and approved by the Zurich Cantonal veterinary office.
Consent for publication
Written informed consent was obtained from all subjects for publication of their individual details and accompanying images in this manuscript. The consent forms are held in the subjects’ notes and are available for review by the EditorinChief.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendices
Appendix 1
Hierarchical Bayesian modelling
The Bayesian inference procedure [30] is summarized as follows. Using the IVIM model
where S _{ n } is a data point measured at the nth bvalue b _{ n } with error term ε _{ n } of Gaussian distribution with variance σ _{ S } ^{2}, the data likelihood reads:
where S = [S _{1}, S _{2}, …S _{ N }]^{T} , g _{ n } = F ∙ exp(−b _{ n } D ^{∗}) + (1 − F) ∙ exp(−b _{ n } D) and N=number of bvalues.
From a Bayesian perspective, nuisance parameters which are of no interest can be marginalized out. Here, the nuisance parameters S _{0} and σ _{ S } ^{2} are integrated out by using a conjugate NormalInverseGamma prior distribution
with g = [g _{1}, g _{2}, …g _{ N }]^{T} and integration over the domain of definition:
The choice of prior allows for an analytic evaluation of Eq. (8). The influence of the prior for the marginalization procedure is diminished by taking the limits of the variance of the normal distribution δ → ∞ and shape and scale parameters of the Inverse Gamma distribution α, β → 0, which encodes a complete lack of prior information. Ignoring proportionality constants which do not depend on the IVIM parameters and are hence not necessary for parameter inference, the marginalized likelihood becomes:
In order to take advantage of the histogram structure of IVIM parameters, a hierarchical prior structure is adopted. First, a multivariate Gaussian is applied to the transformed IVIM parameters. Those transformations are mapping the domain of definition to the field of real numbers:
The multivariate normal distribution considers heterogeneity across the regionofinterest and models correlations between the parameters via the covariance matrix Σ _{ μ }:
The mean across the regionofinterest of the parameter θ _{ i } = [d _{ i }, f _{ i }, d ^{∗} _{ i }] of voxel i is \( \boldsymbol{\mu} =\left[{\mu}_d,{\mu}_f,{\mu}_{d^{\ast }}\right] \). Jeffrey’s prior [31] is used for the hyperparameters μ and Σ _{ μ }:
This prior describes a high probability for a large determinant of the Fisher Information I(Σ _{ μ }) and hence in the considered case of a multivariate normal distribution for a small determinant of the parameter covariance matrix Σ _{ μ } in the regionofinterest. Therefore a “shrinking” of parameter estimates towards the mean of the distribution can typically be observed. The posterior is then given by
The total number of voxels considered here is M. The parameter independent data evidence p(S _{1 : M }), is not required for the inference procedure. Expectation values under the posterior in Eq. (13) are calculated for example as
and analogously for other parameters or quantities.
In order to determine variance under the posterior as uncertainty measure for the LSQ method, a flat prior according to the box constraints in Eq. (4) is chosen instead of the hierarchical prior structure [30].
Markov chain Monte Carlo (MCMC) implementation
The integration in Eq. (14) cannot be performed analytically and therefore a MCMCbased numerical approach is implemented [30]. The expectation value of a parameter (for example d) is approximated by using N _{ s } samples \( {d}_i^{(j)} \)from a Markov chain output:
Appendix 2
Infarcted porcine heart
Both LSQ and BSP approaches were compared in an animal model of acute myocardial infarction. Obtained diffusion weighted data from a single female pig (65 kg, heartrate 68 beats/min) using secondorder motioncompensated diffusion gradients [17] were evaluated. The imaging parameters were: inplane resolution: 2.4 × 2.4 mm^{2}, slice thickness: 10 mm, one apical slice, FOV: 230 × 120 mm^{2}, acquired kspace lines: 49, TR/TE: 2 RR/93 ms, trigger delay 50% peak systole, flip angle: 90°, 8 signal averages and 6 diffusion encoding directions [48]. The 15 optimized diffusion encoding strengths were taken from [46] with a maximum bvalue of 740 s/mm^{2}. The apical myocardial infarct was induced by a permanent distal ligation of the left anterior descending coronary artery (LAD). The animal was anesthetized by a constant dose of Propofol (1.0 ml/kg/min) during surgery and the scan. The experiment was performed in adherence to the Swiss law of animal protection and approved by the Zurich Cantonal veterinary office.
The IVIM parameter maps are shown together with dynamic contrast enhanced (DCE) first pass perfusion images in Fig. 11. The infarcted area can be delineated in the septal area (red arrow). The LSQ and BSP derived diffusion coefficients D are similar across the whole LV. The perfusion fraction F is reduced in the septal area for both methods, with outliers present for LSQ. The pseudodiffusion coefficient D ^{∗} is also clearly reduced in this area for BSP but shows insensible high values for LSQ due to the optimization reaching the fit constraints. If the blood flow related [25, 27] product of F × D ^{∗} is considered, the LSQ derived product does not allow to delineate the infarcted region due to many outliers. In contrast, the septal region in the BSP derived map shows reduced IVIM perfusion parameters.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Spinner, G.R., von Deuster, C., Tezcan, K.C. et al. Bayesian intravoxel incoherent motion parameter mapping in the human heart. J Cardiovasc Magn Reson 19, 85 (2017). https://doi.org/10.1186/s1296801703911
Received:
Accepted:
Published:
Keywords
 Cardiac diffusion imaging
 motioncompensated diffusion weighted spinecho
 intravoxel incoherent motion
 Bayesian inference
 perfusion