An assessment of Bayesian IVIM model fitting
Oscar Gustafsson1,2, Mikael Montelius1, Göran Starck1,2, and Maria Ljungberg1,2

1Department of Radiation Physics, Institute of Clinical Sciences, The Sahlgrenska Academy, University of Gothenburg, Gothenburg, Sweden, 2Department of Medical Physics and Biomedical Engineering, Sahlgrenska University Hospital, Gothenburg, Sweden


Bayesian model fitting has been proposed as an alternative to the commonly used least squares fitting of the IVIM model. In this work we used Monte Carlo simulations to study the convergence of a Markov Chain Monte Carlo implementation of Bayesian model fitting and compared the resulting model parameters to two least squares model fitting methods. We saw that the convergence of the Bayesian model fitting procedure was affected by noise and compartment sizes. Bayesian model fitting was beneficial for the diffusion coefficient and the perfusion fraction, especially at low SNR

Introduction and purpose

Intravoxel incoherent motion (IVIM) has been proposed as a method for monitoring diffusion and perfusion tissue properties1. The biexponential model used for IVIM can be expressed as: $$ S=S_0((1-f)e^{-bD}+fe^{-b(D+D^{*})})$$ where D is the diffusion coefficient, D* is the pseudodiffusion coefficient, f is the perfusion fraction, S0 is the signal without diffusion weighting and b characterizes the diffusion weighting. However, it was early noticed that least squares fitting resulted in noisy parameter estimates2. Bayesian model fitting was then proposed as a more stable method2 and has since been further developed for in vivo applications3. The aim of this work was to study the convergence of Bayesian model fitting parameters and to assess the potential gain in using Bayesian fitting compared to other commonly used methods.

Materials and methods

The Bayesian model fitting was performed using Markov Chain Monte Carlo with a similar set up as Orton et al3. The noise variance was analytically marginalized using a Jeffery prior and assuming a Gaussian noise distribution4. Flat priors within parameter limits were used for all other parameters (0<D<3μm2/ms, 0<D*<1000μm2/ms, 0<f<1, 0<S0). Before the start of the sampling in the parameter probability space, the step length parameters were optimized every 100th iteration for B iterations using the assignment: $$w_{ij}:=w_{ij} \frac{101}{2(101-A_{ij})}$$ where wij is the step length parameter of the i:th model parameter (i=[1,2,3,4]) in the j:th voxel and Aij is the number of accepted samples during the last 100 iterations. The aim of this assignment formula is to produce step length parameters that give samples of which approximately 50% will be accepted3. All step length parameters were initialized to one tenth of the corresponding model parameter's start value. After the last update of the step length parameters another B iterations were performed to reduce the effect of the initial values followed by N iterations with sampling every iteration. To summarize the sampled parameter probability distributions, the mean was calculated for each individual model parameter.

The convergence of the step length parameters as well as the actual model parameters (using converged step length parameters) was studied by monitoring the evolution of the parameters as the number of updates/iterations increased.

The effects of noise levels and compartment sizes were studied by varying the SNR (10, 20, 40) and the perfusion fraction (0.05, 0.15, 0.25). Other parameters used in the simulations were: D = 1μm2/ms, D* = 50μm2/ms, S0 = 100 and b = [1.5,5,10,20,35,50,75,100,200,400,600,800] s/mm2. 1000 voxels were simulated.

The performance of the optimized Bayesian model fitting was compared to two commonly used methods: non-linear least squares fitting (NLLS) and a two-step procedure (2-step) including, (1) a monoexponential fit using b-values above a certain threshold (200 s/mm2 in this study) to get the estimate of D, and (2) another monoexponential fit using all b-values to get the estimates of D* and f, while D is fixed. The 2-step method assumes that the signal from the perfusion compartment is negligible above the specified threshold. 10,000 voxels were simulated. The results were summarized as the absolute error of each parameter p defined as: $$ \sigma_p = \sqrt{\frac{\sum^P_{i=1}{(p_i-p)^2}}{P}} $$ where pi is the fitted parameter value in the i:th voxel, p is the true parameter value and P is the total number of simulated voxels.

Results and discussion

As expected, the step length parameters converged to different values for different levels of SNR and perfusion fraction (figure 1). Convergence was seen for all step length parameters well before 20 updates. 20 updates were therefore used in the subsequent fitting.

Similarly, the model parameters' sample distribution mean was found to be converged before 15000 iterations, which was used in the model fitting method comparison (figure 2).

For most noise levels and perfusion compartment sizes the Bayesian method proved to be superior for the diffusion coefficient and the perfusion fraction (figure 3). However, for D* the 2-step method was the best for almost all cases. The discrepancy compared to previous results for the D* parameter2, where the Bayesian method was superior to NLLS, is most likely due to the use of different prior distributions. When little information about a parameter is contained in the data, the posterior parameter distribution is basically a copy of the prior distribution.


Convergence of Bayesian model fitting was shown to be affected by noise level and compartment size. The best performing fitting method was found to vary as the noise level and compartment sizes were changed. At noise levels commonly seen in practice, Bayesian model fitting provides an attractive alternative without ad hoc thresholds and the possibility of acquiring non-parametric parameter uncertainties.


This research was supported by the Swedish Research Council, the Swedish Cancer Society, the King Gustaf V Jubilee Clinic Cancer Research Foundation in Gothenburg, Sweden and grants from the Swedish state under the LUA/ALF agreement.


1. Le Bihan et al. 1988 Separation of Diffusion and Perfusion in Intravoxel Incoherent Motion MR Imaging. Radiology

2. Neil et al. 1993 On the Use of Bayesian Probability Theory for Analysis of Exponential Decay Data: An Example Taken from Intravoxel Incoherent Motion Experiments. Magnetic Resonance in Medicine

3. Orton et al. 2014 Improved Intravoxel Incoherent Motion Analysis of Diffusion Weighted Imaging by Data Driven Bayesian Modeling. Magnetic Resonance in Medicine

4. Bretthorst et al. 2005 Exponential Parameter Estimation (in NMR) Using Bayesian Probability Theory. Concepts in Magnetic Resonance Part A


Figure 1. Evolution of the median step length parameter as the number of updates increases

Figure 2. Evolution of the median parameter mean value as the number of iterations increases

Figur 3. Absolute error for the three methods used for model fitting

Proc. Intl. Soc. Mag. Reson. Med. 24 (2016)