Our aim was to evaluate different fitting methods for intravoxel incoherent motion imaging model (IVIM) (1) and compare these methods with the monoexponential, kurtosis, and stretched exponential models/functions in the terms of fitting quality, repeatability, and prediction of prostate cancer (PCa) aggressiveness

Eighty-one patients with histologically confirmed PCa underwent two MR examinations on the same day performed using a 3T MR scanner (Ingenuity PET/MR, Philips, Cleveland, USA). The DWI was performed using a single shot SE-EPI sequence, monopolar diffusion gradient scheme, and the following parameters: TR/TE 1394/44 ms, FOV 250x250 mm2, acquisition matrix size 124x124, reconstruction matrix size 256x256, slice thickness 5.0 mm, no intersection gaps, diffusion gradient timing (Δ) 21.204 ms, diffusion gradient duration (δ) 6.600 m, SENSE factor of 2, partial-Fourier acquisition 0.69, SPAIR fat suppression, NSA 2, b values 0, 2, 4, 6, 9, 12, 14, 18, 23, 28, 50, 100, 300, 500 s/mm², acquisition time 3 minutes 45 seconds. The mean signal intensity of squared shaped ROI (4.89x4.89x5.00 mm³), placed in the center of PCa area, peripheral zone (PZ), and central gland (CG), was fitted. The IVIM biexponential equation (Eq. 1) was fitted using the following five different fitting methods:

$$S(b)=S_{0}(fe^{-bD_{p}}+(1-f)e^{-bD_{f}})$$ Eq. 1

1. “Full method”: All four parameters (S₀, D_p, D_f, f) were derived using least square fitting method, in-house written C++ code, utilizing Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm (2) in dlib library (3).

2. “Segmented method”: In the first step, the monoexponential equation (Eq. 2) was used to derive D_f parameter value by fitting signal intensities in the range 100 - 500 s/mm². $$S(b)=S_{0}e^{-bD_{f}}$$ Eq. 2

In the second step, the D_f parameter value from the first step was inserted into the biexponential equation and the remaining three parameters (S₀, D_p, f) were fitted.

3. “Over-segmented method” (4, 5): The first step consisted of fitting signal intensities of b values equal to or higher than 100 s/mm² with the monoexponential model (Eq.2) identically to the first step of the “segmented method”. In the second step, the extrapolated signal of the fitted monoexponential model was used to estimate f according to the equation 3:

$$f=(S_{0}-intercept)/S_{0}$$ Eq.3

, where the intercept is the S₀ estimated from the eq. 2

In the last step the D_f and f parameter values from the first and second step, respectively, were inserted into the biexponential equation (Eq.1) and the remaining 2 parameters (S₀, D_p) were fitted.

4. “Semi-continuous multi-exponential method ” (6-8): This fitting method utilizes Non-negative least Squares (NNLS). The assumption is that the decay curve is composed of multiple mono exponential components each with different fraction of contribution to the decay curve extracted by NNLS. Accordingly, arbitrarily number of coefficient between 0.1 and 1000 μm²/ms were chosen to derive the diffusion distribution spectrum of the signal decay curve. Fraction f was determined from the spectrum by calculating the ratio of the integral between 10-100 μm2/ms and the total integral. $$S(b)=\sum_i^N f_{(ADC_{i})}.e^{-bADC_{i}}$$ Eq.4

5. “Simplified IVIM Model” (9): The biexponential function was modeled using delta function (Dirac delta function, d). Diffusion signal decay reduces to a mono-exponential function for all non-zero b-values: $$S(b)=S_{0}(f\delta(b)+(1-f)e^{-bD_{f}})$$ Eq. 5

In addition to the IVIM fitting methods the following mathematical models/functions were fitted:

1. Monoexponential model (10): $$S(b)=S_{0}(e^{-bADC_{m}})$$ Eq. 6

2. Kurtosis model (11): $$S(b)=S_{0}(e^{-bADC_{k}+\frac{1}{6}b^{2}ADC_{k}^{2}K})$$ Eq. 7

3. Stretched exponential model (12, 13): $$S(b)=S_{0}(e^{-(bADC_{s})^{\alpha}})$$ Eq. 8

The fitting quality was evaluated using corrected Akaike information criteria difference (ΔAICc) (14) while the repeatability of the fitted parameters was evaluated using coefficient of repeatability (CR) and Intraclass Correlation Coefficient (ICC) values (15), specifically ICC(3,1). Receiver operating characteristic curve (ROC) analysis was used to evaluate ability of the fitted parameters (17 parameters in total) to discriminate PCa with Gleason score of 3+3 from those with Gleason score of >3+3. Spearman correlation coefficient (ρ) values were calculated between the fitted parameters and the Gleason score groups (n=3).

Based on ΔAICc the monoexponential model was the preferred model over all of the remaining models/functions and IVIM fitting methods in PCa, PZ and CG (Figure 1). The CR, ICC(3,1), AUC, and ρ values of ADC parameters (ADC_m, ADC_s, ADC_k) were similar to the D_f parameters estimated using all of the IVIM fitting methods (Figure 3). In contrast all f and Dp parameters demonstrated low repeatability (CR, ICC(3,1) values) and diagnostic performance (AUC and ρ values).Monoexponential model demonstrated the highest repeatability and clinical performance in the regions-of-interest based analysis of PCa DWI obtained using b values in the range of 0-500 . s/mm².