An optimized b‐value distribution for triexponential intravoxel incoherent motion (IVIM) in the liver

To find an optimized b‐value distribution for reproducible triexponential intravoxel incoherent motion (IVIM) exams in the liver.


| INTRODUCTION
In 1988, Le Bihan et al 1 introduced the intravoxel incoherent motion (IVIM) model to separate perfusion from diffusion effects within an MRI measurement. The measured signal can be described by a biexponential function: The perfusion fraction f quantifies the relative signal of incoherently flowing blood, the pseudo-diffusion coefficient D * represents blood perfusion, and D is the tissue diffusivity.
compartments because the then expected, B 0 dependency was not observed. 15 A possible explanation is that the IVIM signal decay curve of the perfusion compartments represents ⟨exp (i v)⟩ = ∫ P (v) exp (i v) dv, where = ∫ G (t) tdt is the flow-weighting generated by the applied gradients G (t), that is, essentially the first gradient moment of the diffusion encoding gradients, and v the velocity. At least in the short-time limit, this should hold true, when blood flow directions do not change during the encoding period. The expectation brackets describe the presence of a velocity distribution P (v) whose shape and origin needs to be further explored. It most presumably has venous, arterial, and capillary contributions. Future studies are needed to explore the nature of P (v) and of the respective blood flow correlation times.
To perform such studies, however, a stable fitting of triexponential IVIM curves is necessary. Estimated fit parameters vary widely across different studies that fitted the triexponential IVIM function (see Table 1). An important aspect to consider is the use of different b-values. Varying the distribution of b-values in the protocol affects the estimates of IVIM parameters. 21,22 In a recent study, b-values down to b ≈ 0.1 s/mm 2 were used, which led to increased determined pseudo-diffusion coefficients, which is in good agreement with the finding of threshold dependent IVIM parameter in the biexponential case. 23 Considering the problem of hardly reproducible pseudodiffusion coefficients 22 and echo time dependent perfusion fractions, 24 the aim of this work was to find optimized b-values that reduce the fit uncertainty in all triexponential IVIM parameters, which are robust at different TEs and for the wide range of pseudo-diffusion coefficients found in liver exams. (2) T A B L E 1 Published triexponential IVIM parameters with a field strength of B 0 = 3T

| Optimization routine
The optimization and all simulations were performed in MATLAB (The MathWorks Release R2017b, Natick, MA). The optimization was performed using the triexponential IVIM Equation 1 as signal equation. To perform the optimization, one must know the typical range of values of the triexponential IVIM parameters. Table 1 provides an overview of currently published values. One thing to note is that D * 2 reported in Riexinger et al 15 is roughly by 1 order of magnitude larger than the D * 2 values reported in the other publications, 5,10,11,13,14 which most likely can be attributed to the lower minimal b-value used in Riexinger et al. 15 The optimization performed here is based on the assumption that these large D * 2 values are a good representation of the real values and that consequently very small b-values, say below 1 s/mm 2 , are needed to sample such large D * 2 values. Triexponential IVIM parameters obtained with such very small b-values have so far been reported only for the liver to our knowledge. Consequently the optimization performed here is based on these liver values.

| Ground truth values
To avoid over-fitting to 1 particular set of triexponential IVIM parameters and to account for the likely spread of values in a patient population, the optimization was based on 27 different combinations of triexponential IVIM values. It is known that the perfusion fraction f found in the biexponential IVIM model depends on TE. 24 A similar effect was observed for f 1 and f 2 , which both were reported to decrease by ~40% in the liver if the TE is reduced from 100 ms to 50 ms. 25,26 Therefore, we assume here that f 1 and f 2 scale identically with the TE and the following values were used in the optimization: The factors 1.0, 0.8, and 0.6 shall account for a TE range of ~45-100 ms. 24,26 The absolute values for f 1 and f 2 (i.e., 17.0% and 13.7%) were set to the mean of the values measured in the 10 first volunteers of Riexinger et al. 15 Concerning the pseudo-diffusion coefficients, we assumed that their true values might be larger or smaller than reported in Riexinger et al 15 given their large fit uncertainty. Hence, D * 1 = 0.085mm 2 ∕s ⋅ {0.6, 1.0, 1.4} and D * 2 = 2.55mm 2 ∕s ⋅ {0.6, 1.0, 1.4} was used and all 9 combinations of these values were probed. Taken together with the 3 combinations for the perfusion fractions, this amounted to 27 combinations of triexponential IVIM parameters. The tissue diffusion coefficient was set to 0.0012 mm 2 /s.

| Evaluation of b-value distributions
The evaluation of the quality of a b-value distribution was based on the summed relative error . The individual relative errors are defined with respect to their ground truth values, for example, for f 1 : Here, f 1,i is the f 1 value obtained from signals that were computed according to Equation 2. Before the fit was performed, Gaussian noise was added to the complex signal to simulate a Rician distribution of the magnitude signal. This was performed N times with a signal-to-noise ratio (SNR) of 30, roughly matching typical in vivo SNRs. 6

| Fitting routine
For the signal fitting, a 2-step approach was used. First, a monoexponential function S (b) = S (0) ⋅ exp (−b ⋅ D) was fitted using only b-values above 200s∕mm 2 and starting points S (0) = 0.8, D = 0.001 mm 2 ∕s to estimate the tissue diffusivity D. Fixing this D value, in the second step, a triexponential signal curve (Equation 2) was fitted to the signal. For this purpose, the signal was normalized by the mean of all b = 0 s/mm 2 signals to ensure S (0) = 1 and to reduce the number of fit parameters. The triexponential fit was performed using the MATLAB function lsqcurvefit (starting points: f 1 = f 2 = 10%, D * 1 = 0.5 mm 2 /s, D * 2 = 10 mm 2 /s; lower bound: f 1 = f 2 = 1%, D * 1 = 0.01 mm 2 /s, D * 2 = 0.01 mm 2 /s; upper bound: f 1 = f 2 = 50%, D * 1 = 5 mm 2 /s, D * 2 = 20 mm 2 /s). The perfusion fractions represent the signal arising from blood. Because only a maximum of 100% of the signal may arise from blood, the sum of upper boundaries was equal to 100%. Because we did not want to favor 1 of the fractions, we set their upper limits to 50%.
The upper limits of the pseudo-diffusion coefficients were set at least 8-fold higher than their expected value based on reference. 15
First, an initial distribution consisting of 6 b-values was generated. Then additional optimized b-values were added consecutively. Figure 1 shows a flow chart of the process. (3)

| Step 1: initial distribution
The initial distribution was chosen by generating 1000 random and unique b-value distributions. Distributions whose fit did not converge were considered as suboptimal and were, therefore, discarded. To ensure that the signal normalization could be performed, at least 1 b-value with b = 0 s/mm 2 was specified. N was initially set to 20. In a second step, all distributions whose σ sum was smaller than twice the minimal σ sum of all considered distributions were tested a second time with N = 40. The distribution with the lowest sum was further optimized by replacing 1 b-value (except b = 0 s/mm 2 ) consecutively with the other allowed b-values (with N = 2000). In this step, to speed up the process, only every third b-value of the predefined 64 b-values were tested first. The 4 b-values around the b-value with minimal σ sum were tested and the distribution with the lowest sum was selected as a new optimized distribution. This step was performed for each b-value and repeated 5 times. The resulting b-value distribution was used as initial b-value distribution.
F I G U R E 1 Optimization scheme: first, an initial 6 b-value distribution was generated. Additional b-values were consecutively added until a set of 100 optimized b-values was obtained

| Step 2: adding of optimal b-values
Starting from the initial b-value distribution, distributions with more b-values were found by iteratively adding 1 new b-value. Each allowed b-value was tested once and the new b-value was set to the b-value with minimal sum using N = 2000. This process was repeated until a distribution of 100 b-values was found.

| Step 3: evaluation
To evaluate the time-normalized fit quality � Three b-value distributions were taken from the literature, 6,10,23 which all featured 16 b-values (Table 2) and were compared to the here optimized first 16 b-values. For each of the distributions, σ sum , f 1 , f 2 , D * 1 , and D * 2 were computed for SNR ranging from 5-200 in steps of 5.
For 1 IVIM parameter set ( f 1 = 17%, f 2 = 13.7%, D = 0.0012 mm 2 /s, D * 1 = 0.085 mm 2 /s, D * 2 = 2.55 mm 2 /s) and 5 different SNRs (30,40,50,60,200), 95% confidence intervals of the fitted IVIM parameters were computed. Therefore, the same fitting procedure as described in "Fitting routine" was used, except for the use of the MATLAB fit-option multistart to generate 99 additional starting points. Because confidence intervals cannot be computed by the use of the MATLAB function multistart, a single fit identical to the multistart fit using only the best start parameters found by the multistart approach, was used. For each IVIM parameter and each number of b-values, the median of the 95% confidence limits estimated by the fit were calculated over all N. Once, the optimized b-values were used (with 6 to 100 b-values). Once, b lit2 was used, which is a b-value distribution taken from indicates a better improvement of fit stability than obtained for linear fits, with the number of acquired b-values n literature featuring 16 b-values. As above, a range of repetitions was used in this case (i.e., amounting effectively to 32, 48, 64, etc. b-values). N = 500 was chosen. Additionally, the Akaike's information criterion (AIC) was computed with 3 SNRs (30, 60, 200) for the optimized b-value distribution (again with N = 500) to compare the bi-and triexponential fit.

| In vivo IVIM measurements
The liver of 3 healthy volunteers were examined on a 3T scanner (Magnetom Prisma; Siemens Healthcare Erlangen, Germany) with an 18-channel body coil and 4 different distributions of 16 Table 2) as well as with the full distribution of | 2101 RIEXINGER Et al.
100 optimized b-values used as reference. The study protocol was approved by the local institutional ethics committee and written informed consent was obtained from all participants.
Diffusion-weighted images were acquired with respiratory triggering, an isotropic voxel size of 4 × 4 × 4 mm 3 , a field of view of 400 × 400 mm 2 , 6 transversal slices, slice gap Effective b-values were subsequently calculated taking into account the impact of imaging gradients (the impact of imaging gradients on the IVIM parameter estimation is shown in Supporting Information Figure S1). This effect is most pronounced for the small b-values. The effective b-values were 0.15, 0.216, 0.298 s/mm 2 for the nominal values b = 0, 0.2, 0.3 s/mm 2 , respectively. In each slice, a region of interest (ROI) was positioned in the first unweighted image, avoiding major blood vessels and only including F I G U R E 5 Median of upper and lower limit of relative 95% confidence intervals (normalized by the true value, that is, for example, for the upper limit of f 1 : (median(upper limit of CI f1 ) − f 1,true )∕f 1,true is plotted) of the optimized b-value set (lines) and multiples of b lit2 (x). The confidence interval of D * 2 estimated with b lit2 were orders of magnitudes higher and therefrom omitted to achieve a better readability | 2103 RIEXINGER Et al.
liver tissue. The signal was calculated using the median intensity over all diffusion directions inside the ROI. It was then normalized to the mean of the unweighted signals. The entire data set, acquired with all 100 optimized b-values, and the 4 individual data sets consisting of 16 b-values each, were fitted with the above described fitting routine (see section "Fitting routine"), which was identical except for the use of the MATLAB fit-option multistart to generate 999 additional starting points and an additional biexponential fit to calculate AIC values and compare the bi-and triexponential fit. Because a reliable ground truth was not available for in vivo measurements, the fitted parameters of the distribution of 100 optimized b-values were used as reference. For each IVIM parameter obtained with the distributions consisting of 16 b-values, the root mean square deviation to the reference values were calculated. Afterward, they were summed to compute a total relative error. Additionally, a pixel-wise evaluation for all slices for 1 volunteer was performed using the same fitting routine. The SNR of each volunteers' liver parenchyma at b = 0 s/mm 2 was estimated by the differential method 27,28 applicable for images acquired with parallel imaging. Two circular ROIs were drawn on each slice into the right liver lobe on the 2 unweighted acquisitions closest in time. The SNR for each ROI was calculated as follows: where S 1 and S 2 are the signal values in each ROI for both time points. The estimated mean SNR of all ROIs in 1 volunteer was defined as the volunteers' SNR. Table 2 and illustrated in Figure 2A This indicates that one gains more than the usual square root of time increase in fit stability found for linear fits; and that one should use more than 6 b-values. A distribution of 100 optimized b-values performed better than repeated acquisitions with distributions of either 6 or 16 optimized b-values. Using the 6 b-value distribution is not recommendable because averaging the data multiple times does decrease the time-normalized fit quality.  Figure 3 shows the comparison of b opt to the b-value distributions taken from the literature. The diffusion coefficient is well-determined by each of the 4 considered 16-b-value distributions, whereas large differences were found when fitting IVIM parameters. b lit2 performed better than b lit1 and b lit3 for estimating f 1 and f 2 , but not as good as b opt . This behavior can be explained by the b-values included in the distribution. The smallest b-value used by b lit1 and b lit3 is 10 s/mm 2 , which seems to be too high for determining the perfusion fractions correctly. For D * 1 , b lit2 and b opt performed equally well and much better than the other 2 distributions. A b-value of 5 s/mm 2 , therefore, appears to be sufficiently small to ensure decent fit precision for D * 1 . For D * 2 , only b opt performed well. Because the initial signal drop, which is linked to the pseudo-diffusion coefficient of the fast compartment, occurs for b < 1 s/mm 2 , b lit1 , b lit2 , and b lit3 were not suited to fit D * 2 . Figure 4 shows the difference of AIC values obtained with bi-and triexponential fit in dependence of the number of b-values taken from the optimized distribution. Negative values in Figure 4 indicate that the triexponential fit represents the data better. If the number of b-values is small, the biexponential fit represents the data better. If more b-values are used, the opposite holds true. With increasing SNR, less b-values are needed to favor the triexponential fit.

The optimized distribution of b-values is listed in
In Figure 5, the medians of the upper and lower limit of the relative confidence intervals of the triexponential fit are shown. Similar to the time-normalized fit quality, the confidence interval decreases with increased SNR and number of b-values. Even for low SNRs, the perfusion fractions f 1 and f 2 , as well as the pseudo-diffusion coefficient D * 1 can be roughly estimated. However, for estimating the fast pseudo-diffusion coefficient D * coefficient is more or less constant and is strongly reduced at 19 b-values. This can be explained by the addition of a high b-value to the distribution.

| In vivo measurement
The estimated SNR was 65, 68, and 72 for the 3 volunteers. Sample image data of 1 volunteer and the used ROI are shown in Figure 6A. The respective signal attenuation curves are plotted in Figure 6B, with the performed bi-and triexponential fit and their corresponding AIC-values. AIC tri − AIC bi is negative, indicating that the triexponential IVIM function is better suited to describe the data than the biexponential function.
The group averaged estimated triexponential IVIM parameters of the optimized data set were: Standard deviations and mean values were calculated among all volunteers and slices. The relative errors and the summed relative error for each investigated b-value distribution, computed with respect to the distribution of 100 optimized b-values, is illustrated in Figure 7. The diffusion coefficient is well determined by all distributions, but the optimized distribution is slightly more inaccurate presumably because of the lower number of high b-values. The perfusion fractions f 1 and f 2 and the fast pseudo-diffusion coefficient D * 2 were determined best by the optimized distribution. The perfusion fractions and D * 1 had the lowest relative errors with b opt , whereas the slow pseudo-diffusion coefficient D * 1 is determined best with b lit2 , which had been optimized for biexponential IVIM. By far, D * 2 (5) D = 0.0011 ± 0.00024mm 2 ∕s, f 1 = 13.1 ± 4.5%, f 2 = 10.8 ± 1.6%, D * 1 = 0.016 ± 0.007mm 2 ∕s, D * 2 = 0.50 ± 0.224mm 2 ∕s. can be determined best with b opt and it dominates, moreover, the summed relative error. The low b-values contained in b opt appear to be crucial for reliably determining D * 2 . IVIM parameter maps could sucessfully be calculated for all slices using the set of 100 optimized b-values. Maps of 1 representative slice are shown in Figure 8 with the corresponding unweighted image. Within the liver, D increases in the left liver lobe, because of the pulsation artifact. 29,30 f 1 seems homogeneously distributed across the whole liver. Bright spots are located in the vicinity of large vessels. f 2 , D * 1 , and D * 2 also exhibit bright spots close to vessels, which are more prevalent than in the f 1 map. Within those spots, the pseudo-diffusion coefficients increase drastically compared to the surrounding tissue. For the slow pseudo-diffusion coefficient, this increase is from 0.006 mm 2

| DISCUSSION AND CONCLUSIONS
The purpose of this study was to find an optimized b-value distribution for triexponential IVIM in the liver. For the biexponential IVIM model, b-value optimizations were already performed for different organs and purposes to increase IVIM parameter reproducibility. 6,[31][32][33][34] Besides the optimization of the b-value sampling, lots of other optimizations for the biexponential IVIM model in the liver, for example, concerning breathing schemes or minimizing the acquisition time were performed. [34][35][36][37][38][39] However, to our knowledge, no comparable optimization has been performed for triexponential IVIM yet. Within this study, we were able to find an optimized distribution of b-values for triexponential IVIM in the liver and demonstrated the advantage of optimized b-values compared to 3 representative b-value distributions used in literature by means of simulations and volunteer measurements. In contrast to b-value distributions used in literature, 6,10,13 our results would suggest that for accurate triexponential IVIM measurements, frequent acquisition of very small b-values (<6 s/mm 2 ) is required. This recommendation is supported by the work of Cohen et al. 21 They showed that the pseudo-diffusion coefficient in the biexponential IVIM model is underestimated, if too few small b-values are acquired, which confirms our finding for the triexponential model. Additionally, pseudodiffusion coefficients exhibit low reproducibility. 40 Therefore, by spending more time on measuring very small b-values, where the signal is highly sensitive to the pseudo-diffusion effects, one can increase the fit precision of triexponential IVIM parameters.
Although the summed relative error within the simulation is dominated by the relative error of the fast pseudo-diffusion coefficient, the other triexponential IVIM parameters were also estimated best by the optimized b-value distribution, except for the slow pseudo-diffusion coefficient, which is estimated best by the b-value distribution optimized for biexponential IVIM by Lemke et al. 6 In our work, only 16 b-values were used for the in vivo measurement for comparison and evaluation of the full optimized b-value distribution. Using more b-values, one can further increase the fit precision by decreasing the confidence interval as shown above. However, the time-normalized fit quality predicts the largest improvement in the range of 10 to 20 b-values for all considered SNRs. Therefore, we would recommend acquiring at least 16 b-values for triexponential IVIM exams of the liver.
The estimated mean pseudo-diffusion coeffiecients D * 1 and D * 2 in the volunteers study were smaller than the ones used as ground-truth values within the optimization. This might be explained by the different slice orientation (axial here, sagittal in Riexinger et al). 15 Another explanation might be the presence of different inflow effects because of the use of respiratory triggering compared to free breathing in Riexinger et al. 15 Free breathing might also lead to inclusion of contributions from larger vessels in the ROI, which presumably leads to increased pseudo-diffusion coefficients. This effect is highlighted in the estimated IVIM maps (Figure 8), where an inhomogenious behaviour of IVIM parameter across the liver and bright spots are clearly observable. These spots are most presumably to be attributed to large vessels. If a ROI-wise evaluation is performed, the amount of larger vessel will have a critical impact on measured pseudodiffusion coefficients. Optimized b-values for the IVIM parameters determined in this study are provided as supplemental material (Supporting Information Table S1). The origin of the observed discrepancy between estimated and groundtruth values is most likely a mixture of issues mentioned above. This illustrates the difficulty of producing triexponential IVIM parameter across different studies and should be the subject of further investigations.
We acknowledge several limitations of this study. First, we did not evaluate the optimized distribution within a clinical study. Although the optimization considered varying IVIM parameters, the optimal b-value distribution for assessment of diseased tissues might differ. Second, the optimized b-value distribution added 1 b-value consecutively instead of optimizing b-value distributions with a fixed number of b-values. Potentially a full optimization carried out for each particular number of b-values might have yielded further improved distributions. Third, in vivo examinations were performed with a comparatively short TE (TE = 45 ms), and

SUPPORTING INFORMATION
Additional Supporting Information may be found online in the Supporting Information section.