Intravoxel incoherent motion (IVIM) has been proposed as a method for monitoring diffusion and perfusion tissue properties¹. 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, S₀ 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 estimates². Bayesian model fitting was then proposed as a more stable method² and has since been further developed for in vivo applications³. 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.

The Bayesian model fitting was performed using Markov Chain Monte Carlo with a similar set up as Orton et al³. The noise variance was analytically marginalized using a Jeffery prior and assuming a Gaussian noise distribution⁴. Flat priors within parameter limits were used for all other parameters (0<D<3μm²/ms, 0<D*<1000μm²/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 w_ij is the step length parameter of the i:th model parameter (i=[1,2,3,4]) in the j:th voxel and A_ij 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 accepted³. 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μm²/ms, D* = 50μm²/ms, S₀ = 100 and b = [1.5,5,10,20,35,50,75,100,200,400,600,800] s/mm². 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/mm² 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 p_i 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.

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* parameter², 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.