The intravoxel incoherent motion model (IVIM) is of interest as it gives a more complete characterization of DWI signal attenuation in a number of normal organs and pathologies [1]. In practice, estimates of the pseudo-diffusion coefficient $$$\small{}D^*$$$ are noisy, so a number of advanced techniques have been proposed: using a fixed value for $$$\small{}D^*$$$ [2], a hierarchical Bayesian model incorporating a histogram model [3], and Markov random field (MRF) models [4]. MRF methods augment the data cost function with a spatial smoothness prior that encourages neighboring voxels to have similar values. MRFs use weighting factors that are typically fixed to subjectively pre-determined values [4,5]. In this abstract, a Bayesian model and associated Markov chain Monte-Carlo (MCMC) algorithm are presented that enable data-driven estimation of the MRF smoothness weights.

To explore the feasibility of the proposed MRF model and MCMC algorithm, and to show that the convergence properties of the algorithm shed light on which parameters of the IVIM model can be reliably estimated from data similar to that obtained in a clinical setting.

MRF model: The data likelihood assumes Gaussian errors between the log-data and the log-transformed IVIM model:$$$\small{}\;\log{}S_i=\log{}S_0+\log\left(\,f\mathrm{e}^{-b_iD^*}+(1-f)\mathrm{e}^{-b_iD}\right)+\eta_i$$$, where $$$\small{}f$$$ is the pseudo-diffusion volume fraction,$$$\,\small{}D^*$$$ and$$$\,\small{}D$$$ are fast and slow diffusion coefficients and $$$\small{}\mathrm{var}(\eta_i)=\sigma_S^2$$$ over the whole ROI/image. The smoothness prior is defined by a product of nearest-neighbour clique terms, each of which is modeled with a Laplace distribution$$$\,\small{}p\!\left(x_i,\,x_j\right)=\frac{1}{2}w_x\exp\left(-w_x\left|x_i-x_j\right|\right)$$$, and the scale term $$$\small{}w_x$$$ is the smoothing weight for $$$\small{}x\in\{f,\,D,\,D^*\}$$$. This distribution is the probabilistic equivalent to the robust L₁-norm [4] which has desirable edge-preserving properties. The posterior distribution is the product of the data likelihood and the smoothness priors for$$$\,\small{}f,\,D\,$$$and$$$\,\small{}D^*$$$. The key innovation introduced here is to estimate the smoothing weights in addition to the voxel IVIM parameters using an MCMC algorithm.

Data and evaluation: Abdominal DWI data were acquired coronally in seven patients with liver tumours on a 1.5T MAGNETOM Avanto (Siemens Healthcare, Erlangen, Germany) using a prototype single-shot EPI sequence and 20$$$\!\small{}\times$$$5mm$$$\;$$$slices, TR/TE=5000/60ms, matrix=128², FOV=400mm, 5$$$\;$$$NSA, 3$$$\;$$$orthogonal directions, phase-partial-Fourier$$$\;$$$7/8, GRAPPA$$$\;$$$2, b-values=0,20,40,60,120,240,480,900s/mm². Individually acquired images were registered [6], and averaged per b-value. A single slice was selected covering the liver, and including the kidneys where visible, and a rectangular portion of the images was cropped to remove the background, see figure 1. Initial values were obtained using least-squares, and the MRF/MCMC algorithm used to generate 10,000 posterior samples.

The principal result is that in all seven patients, estimates of $$$\small{}w_{D^*}$$$ did not converge, but instead diverged to large values, which leads to unrealistically smooth parameter maps. Therefore, for the remainder of the evaluation,$$$\,\small{}w_{D^*}$$$ was fixed to 5. This was chosen to be at the lower end of the estimates of $$$\small{}w_D$$$ and$$$\,\small{}w_f$$$ obtained in figure 2, in order to avoid over-smoothing the $$$\small{}D^*$$$ maps. For the Laplace clique term with weight$$$\,\small{}w$$$, the average percentage change between neighbouring voxels is$$$\,\small{}100\%\times(w-1)^{-1}$$$, which is around 25% for$$$\,\small{}w_{D^*}=5$$$. The example in figure 1 shows the least-squares fit (used as initialisation), parameter maps from a single MCMC sample (showing the level of image smoothness determined by the estimated MRF weights), the mean of the post-burn-in MCMC samples (showing the final MRF estimates), and the coefficient of variation of the MCMC samples (computed from the standard deviation of the MCMC samples). With$$$\small{}\,w_{D^*}\,$$$fixed,$$$\small{}\,w_f$$$ failed to converge in one patient, and$$$\small{}\,w_D\,$$$failed to converge in a different patient. Posterior distributions for the smoothing weights for all cases which converged are shown in figure 2, and the mean$$$\small{}\;(\pm$$$sd$$$)\,$$$were$$$\small{}\,w_D=14.9\pm2.6\,$$$and$$$\small{}\,w_f=8.4\pm2.0$$$.

For the Bayesian model given above, estimation of the smoothing weights is essentially driven by a balance between spatial heterogeneity and parameter uncertainty. For these data, the convergence of$$$\small{}\,w_D\,$$$and$$$\small{}\,w_f$$$ in almost all patients indicates that the uncertainty (given the image noise and b-value support) in$$$\small{}\,D\,$$$and$$$\small{}\,f$$$ is sufficiently low that the degree of spatial heterogeity can be determined from the data. However, the divergence of$$$\small{}\,w_{D^*}$$$ leads to very smooth $$$\small{}D^*$$$ maps, which is equivalent to a single value of$$$\small{}\,D^*$$$ for the whole image. This is consistent with the known difficulty in estimating$$$\small{}\,D^*$$$, and is driven by poor data support for the pseudo-diffusion coefficient of the IVIM model. Despite this finding, it is clear that a single value for$$$\small{}\,D^*$$$ is not appropriate when modelling complex anatomy, which motivates the use of a fixed value for$$$\small{}\,w_{D^*}$$$. To account for image acquisition differences (in particular the voxel dimensions), it would be preferable to estimate the smoothing weights independently for every patient, but where this is impractical, weights can be estimated from a subset of patients and the average used for subsequent patients.