An assessment of Bayesian IVIM model fitting

Oscar Gustafsson^{1,2}, Mikael Montelius^{1}, GĂ¶ran Starck^{1,2}, and Maria Ljungberg^{1,2}

The Bayesian model fitting was performed using Markov
Chain Monte Carlo with a similar set up as Orton et al^{3}. The noise variance was
analytically marginalized using a Jeffery prior and assuming a Gaussian noise distribution^{4}. Flat priors within parameter limits were used for all other
parameters (0<D<3μm^{2}/ms, 0<D*<1000μm^{2}/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^{3}. 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^{2}/ms, D* =
50μm^{2}/ms, S_{0} = 100 and b = [1.5,5,10,20,35,50,75,100,200,400,600,800] s/mm^{2}. 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

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^{2}, 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.

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
Proc. Intl. Soc. Mag. Reson. Med. 24 (2016)

2080