Influence of saturation effects on biexponential liver intravoxel incoherent motion

Studies on intravoxel incoherent motion (IVIM) imaging in the liver have been carried out with different acquisition protocols. The number of acquired slices and the distances between slices can influence IVIM measurements due to saturation effects, but these effects have often been disregarded. This study investigated differences in biexponential IVIM parameters between two slice settings.


INTRODUCTION
Intravoxel incoherent motion (IVIM) is a concept in DWI that considers signal contributions stemming from both the tissue and blood. 1 The so-called IVIM effect in the liver becomes noticeable as a strong signal decay for small b values (roughly b < 150 s∕mm 2 ). Molecular diffusion in tissue at b values smaller than approximately D −1 is usually approximated well by a monoexponential signal decay, where S 0 is the signal without diffusion weighting, b is the diffusion weighting, and D is the tissue diffusion coefficient.
Including the perfusion compartment, the signal attenuation is commonly described by a biexponential function as follows 1,2 : Here, the blood perfusion is characterized by the perfusion fraction f and the pseudodiffusion coefficient D * . The difference in magnitude between D and D* (D* is usually one magnitude larger than D) enables differentiation between the blood and tissue compartments. A variety of acquisition parameters have been used in liver IVIM studies. For instance, a strong dependence of f on the TE was observed. 3,4 Similar to other imaging parameters, the slice distance used and the number of slices have varied considerably among studies. [3][4][5][6][7][8][9][10] For instance, Riexinger et al. 5 used four slices with a spacing of 4 mm between slices, whereas Cercueil et al. 8 performed image acquisition with 20 slices and a small slice distance of 1.4 mm. Lemke et al. 4 used 14 slices and even smaller gaps between slices (0.25 mm). These differences raise the question of whether the obtained IVIM parameters are comparable among studies or whether correction factors are necessary, similar as for the TE. 11 Saturation effects on flowing magnetization may be relevant in this regard. Such saturation effects are widely used in MRI angiography, but their potential presence has been mostly neglected in IVIM studies. Saturated blood magnetization might travel from one slice into adjacent, subsequently measured slices. Consequently, f may decrease.
The aim of this study was to reveal whether such saturation effects are of importance in IVIM examinations of the liver. To this end, IVIM parameters were compared between a setting with a high number of slices and small slice gaps and a setting with a low number of slices and large slice gaps.

Volunteers
Diffusion-weighted images of the liver of 15 healthy subjects with no known history of liver diseases were acquired. The volunteers were scanned in a supine position. The age of the volunteers was between 21 and 30 years (mean age: 25.8 years); five were female and 10 were male. The study was approved by the local ethics committee. All volunteers provided a written declaration of informed consent.

Data acquisition
All measurements were performed using a 3 T scanner (Magnetom Prisma, Siemens Healthineers, Erlangen, Germany) with a single-refocused diffusion-weighted echoplanar imaging product sequence. In this study, the sequence was adapted with the Siemens sequence programming software IDEA(VE11E) (Siemens Healthineers, Erlangen, Germany), allowing for a b-value step size of 1 s/mm 2 (instead of the default step size of 50 s/mm 2 ). An 18-channel body coil and an integrated 32-channel spine coil were used for signal detection.
The b values were optimized based on the "b high " b-value set by Lemke et al. 12 To keep the acquisition time under 20 min, the total number of b values was limited to 16. All b > 800 s∕mm 2 were lowered to 800 s/mm 2 to avoid kurtosis effects. For the final evaluation, b values were corrected by calculating them numerically for each diffusion direction using the POET Siemens sequence simulation tool (Siemens Healthineers, Erlangen, Germany) and taking into account all imaging gradient pulses.  12 and those used in our study is given in Table S1. For each b value, except for b = 0, data were acquired in three orthogonal diffusion directions, (−1, 0, 0), (0, 0, −1), and (0, 1, 0); the directions are stated in the scanner coordinate system.
All measurements were performed in free breathing mode. TE was set to 60 ms and TR to 3600 ms. This

F I G U R E 1
Representative images and ROIs (red) of one slice at different b values. Note the slight shift of the liver position due to respiratory motion. ROIs were adapted accordingly. Imagesare windowed equally. ROI, region of interest.
TR represented the lowest possible value for the setting with a high number of slices. The other acquisition parameters were the slice thickness (5 mm), in-plane voxel size (4 × 4 mm 2 ; interpolated to 2 × 2 mm 2 ), FOV (400 × 400 mm 2 ), acquisition bandwidth (2942 Hz/Px), echo spacing (0.4 ms), and phase partial Fourier factor (6/8). The following vendor-provided data correction algorithms were used: "dynamic field correction" correcting for eddy current-induced distortions, "2D distortion correction" correcting image distortions arising from gradient nonlinearities, "pre-scan normalize" compensating for surface coil flare, and "raw filter" of "medium" strength to reduce Gibbs ringing. The parallel imaging method GRAPPA with an acceleration factor of 2 and 24 reference lines was used. Spectral-attenuated inversion recovery was used for fat suppression. All images were acquired in transversal orientation.
Two different acquisition settings were used. First, data were acquired in a minimal saturation setting. For this, four slices were placed over the entire liver with a spacing that was large enough to reduce saturation effects to a large extent. A distance of 25 mm between slices was used. This slice setting is referred to as the "few slices setting" (FS).
Ideally, a single slice would have been measured for each subject to eliminate any saturation effect. However, this would have impacted the quality of the statistical analysis.
In the second setup, referred to as the "many slices setting" (MS), the slice gap was 1 mm, and the entire liver was covered for a total of 24 to 27 slices depending on the size of the liver.
The total acquisition time for both slice settings added up to approximately 15 min.

Quality check
An initial quality check was performed before processing the data further. If the overall image quality was deemed insufficient upon visual examination, all data of the specific slice was discarded. If an image was affected by pulsation artifacts [13][14][15] expanding into the right liver lobe, the regions of interest (ROIs) of all images were reduced to avoid including the region that was visually affected by the artifact.

Data evaluation
For each slice, a ROI was drawn on a b = 1 s/mm 2 image in the liver parenchyma of the right liver lobe sparing major vessels and was copied to all other images ( Figure 1). This was performed using the Medical Imaging Interaction Toolkit v2021.10 (German Cancer Research Center, Heidelberg, Germany). Due to the respiratory motion of the liver, the ROIs were adapted if necessary to solely cover the liver tissue for all b values. The left liver lobe was deliberately excluded, as it is more prone to pulsation artifacts. [13][14][15] The median value of all voxel signals inside the ROI was then computed separately for each slice, b value, and diffusion direction using Python 3.8. The median was used instead of the mean because it is more robust against outliers. The median signal value of each slice was then normalized to the respective signal value at b = 0 and averaged arithmetically over all slices for each b value and diffusion direction. The signal values were then averaged over the three diffusion directions.
To obtain the IVIM parameters, fitting was performed as a segmented fit. 3,[16][17][18][19][20][21] First, the tissue diffusivity D and the intercept a of the monoexponential signal decay with the ordinate were determined using signals acquired with b values greater or equal to 150 s/mm 2 . A monoexponential signal equation and the Levenberg-Marquardt algorithm were used. Signal values at 800 s/mm 2 were acquired twice and were thus double-weighted in the fit. Starting values for the fit for a and D were the maximum signal value (considering all b values) and 1 ⋅ 10 −3 mm 2 ∕s, respectively.
In the second step, the estimation of f was performed, albeit in a different fashion than previously described 3,16-21 because the data acquired for this study suggested a triexponential IVIM curve for low b values (Figure 2A,B). For this reason, the biexponential IVIM equation did not represent the data very well, and f could not be determined with the commonly used biexponential fit approach. 3,[16][17][18][19][20][21] To obtain the perfusion fraction, f is simply determined by 1 − a. D * and S 0 were determined using a biexponential fit while treating D and f as fixed parameters. Here, bounds were provided for D * to avoid physically non-meaningful parameters while still leaving enough margin to not appreciably affect the fit. These bounds restricted D * to between 0 and 50 ⋅ 10 −2 mm 2 ∕s and S 0 to between 0 and 1.4. An alternative approach would be to fix S 0 for the biexponential fitting. However, we chose the current method as it is Representative biexponential, triexponential IVIM and monoexponential fit. IVIM curve for all b values (A) and enlarged view for small b values (B). At low b-values, a triexponential decay describes the data better. Data points for which the biexponential fit curve slightly deviates are marked with blue errors. The dashed lines indicate the axis at b = 0 and the cutoff b value at 150 s/mm 2 separating the bi-from the monoexponential region for fitting D (C) . The intersection of the monoexponential fit with the b = 0 axis (depicted in green) is used for determining the perfusion fraction without utilizing a biexponential fit. less prone to outliers for D * (see Appendix S1), and it also is in accordance to Führes et al. 3 where biexponential fitting with S 0 has been applied. Additionally, a triexponential fit was performed for visualization purposes (Figure 2) as described by Führes et al. 3

Statistics
The Shapiro-Wilk test was performed to test for normal distribution of the IVIM parameters. For normally distributed data, Student's t test for paired samples was applied. For non-normally distributed data, the Wilcoxon signed-rank test was used. The significance level for the statistical tests was set to 0.05. To investigate linear correlations between datasets of the two slice settings, the Pearson correlation coefficient was computed. All statistical analyses were performed in Python 3.8.

RESULTS
No dataset failed the initial visual quality check. The ROIs in the right liver lobe were reduced in size in approximately 50% of all datasets due to pulsation artifacts extending into the right liver lobe. In Figure 1, representative images are shown for one slice of one volunteer at different b values. The respective ROIs are also shown. Due to the respiratory motion of the liver, slight differences can be observed among the images.
A representative signal curve with an extrapolated monoexponential fit and the intercept with the b = 0 axis is depicted in Figure 2C. Figure 3A shows the boxplots of the biexponential IVIM parameters for both slice settings. The mean values, SDs, median values, first (Q1) and third quartiles (Q3), and results of the statistical analysis are listed in Table 1. No noteworthy difference between the slice settings was observed for D and D * . For f , the mean and the first and third quartiles were reduced for the MS setting, but the median was not.
The data points of D were not normally distributed and thus the Wilcoxon signed-rank test was applied. For f and D * , Student's t test for paired samples was used, as the Shapiro-Wilk test indicated normally distributed data in both slice settings. With P values of 0.934, 0.177, and 0.977, respectively, D, f , and D * were not significantly different between the FS and MS settings. Figure 3B presents the scatter plots for D, f , and D * . The Pearson correlation coefficients for the IVIM parameters were 0.813 (D), 0.478 (f ), and − 0.022 (D * ), indicating a moderate to strong positive correlation between the slice settings for D and f and a negligible correlation for D * . 22,23

DISCUSSION
In this work, the biexponential IVIM parameters of two slice settings were compared to examine the influence of T A B L E 1 Mean, SD, median, first quartile, and third quartile of the biexponential IVIM parameters for the two slice settings, the applied statistical tests, and the respective P values Note: D few and D many are presented as 10 −3 mm 2 ∕s; f few and f many are presented as %; D * few and D * many are presented as 10 −2 mm 2 ∕s. Abbreviations: D, tissue diffusion coefficient; D*, pseudodiffusion coefficient; f , perfusion fraction; IVIM, intravoxel incoherent motion. saturation effects. No significant changes were observed for the IVIM parameters in the two settings.
A meta-analysis by Li et al. 24 that included data from 27 healthy liver studies was used to compare the values of our IVIM parameters with previously published values. The overall median of D was (1.09 ± 0.17) ⋅ 10 −3 mm 2 ∕s. This is roughly in line with our median values of 1.14 ⋅ 10 −3 mm 2 ∕s (FS) and 1.20 ⋅ 10 −3 mm 2 ∕s (MS). The median of f was (22.4 ± 8.5)% , which is lower than the values measured in our study ( 28.1% (FS) and 28.6% (MS)). We attribute this difference to the fitting approach of f ; in our study, we used an extrapolation approach different from the commonly used biexponential fitting procedure.
The reported median value of D * is (7.06 ± 3.10) ⋅ 10 −2 mm 2 ∕s . 24 This is slightly lower than the values found in this study ( 8.43 ⋅ 10 −2 mm 2 ∕s and 8.33 ⋅ 10 −2 mm 2 ∕s ). However, taking the large measurement uncertainties into account, ours are relatively close to those found in the previous studies. The values of D * varied widely among previous studies. 5,9,25,26 For example, Riexinger et al. 5 reported a median value of (26.0 ± 2.6) ⋅ 10 −2 mm 2 ∕s at 3 T. Lee et al. 9 performed two volunteer examinations of the liver. The mean values of D * in the right liver lobe were (13.58 ± 7.82) ⋅ 10 −2 mm 2 ∕s in the first examination and (13.30 ± 8.15) ⋅ 10 −2 mm 2 ∕s in the second examination. Jerome et al. 25 reported a much lower mean value for the liver, at (4.4 ± 0.5) ⋅ 10 −2 mm 2 ∕s . The results of Leporq et al. 26 lie in between, at (8.02 ± 2.33) ⋅ 10 −2 mm 2 ∕s .
In general, the values of D and D * found in the current study are in line with those of previous studies, whereas f was larger than in previous studies.
The fact that D did not show a significant dependency on the slice setting was expected because the slice setting likely does not alter the composition or properties of the liver tissue. The two slice settings both covered the entire liver; thus, variations in the physiological properties of different liver regions should have a similar impact on D.
Similarly, there was no significant difference in D * between the two settings. Large measurement uncertainties occurred in both settings, for example, the interquartile ranges were 7.34 ⋅ 10 −2 mm 2 ∕s (FS) and 4.00 ⋅ 10 −2 mm 2 ∕s (MS), which gives a relative deviation of 87% (FS) and 48% (MS) with regard to the median values. These results indicate that the potential influence of the slice setting on the pseudodiffusion coefficient was probably not measurable in this study due to the large fluctuation of values. A potential source for such a difference could be that faster flowing blood could be saturated more often because it travels more efficiently to adjacent slices. The pseudodiffusion coefficient describes both faster and slower flowing blood. If faster flowing blood is more saturated, then D * might become smaller. Such an effect was not observed, however.
In addition, f did not change significantly between the two slice settings. The mean f was slightly higher for FS (29.7%) than for MS (27.7%), and also Q1 and Q3 slightly decreased, indicating the potential presence of a weak saturation effect, whereas the median values were nearly equal. Student's t test confirmed this result by indicating no significant difference. This proves that saturation effects do not noticeably impact f . In Appendix A, we outline a biophysical model to interpret these results ( Figure 4).
The saturation effects may have a larger influence in the MS setting when much shorter TRs are used. A lower TR results in a reduced blood signal due to the relatively long T 1 time of blood [27][28][29][30][31][32] in comparison to liver tissue. On the other hand, blood could only flow into a few adjacent slices at shorter TR, which could reduce the saturation effect. These circumstances must be validated experimentally in future studies.
In previous studies, the following TR were used, for example: 1500 ms, 33 2000 ms, 26 2500 ms, 5 2600 ms, 7 4000 ms, 25 5000 ms, 34,35 6800 ms, 36 and one respiratory cycle. 3,8 These TR varied between 1500 and 6800 ms (and could exceed 7000 ms depending on the duration of one respiratory cycle). The majority of studies were performed with roughly equal or longer TRs compared with our

F I G U R E 4
Boxplots of the biexponential IVIM parameters for the FS and MS settings. Blue circles each represent the IVIM parameters of one volunteer. Outliers are additionally denoted as red "+"-signs. The red lines indicate the respective medians and dark red "x"-signs show the respective mean values. Boxes mark the IQR, and whiskers indicate all data lying within 1.5 × IQR. FS, few slices; IQR, interquartile range; MS, many slices. study (TR = 3600 ms). However, two studies stand out: Luciani et al. 33 and Leporq et al. 26 Our findings indicate that the variance in the TRs used does not prevent a straightforward comparison of the quantitative IVIM parameters obtained in these studies. Because the influence of TR on the quantitative liver IVIM parameters appears to be quite subtle, a direct comparison of the quantitative values found in these studies with those of our study would be unable to clarify the size of this influence. Other parameters differ as well (e.g., b values), whose influence most presumably outweighs any TR dependency effect. Using the model outlined in Appendix A, we estimate that the saturation effects might indeed be relevant for decreased TR. The increasing availability of simultaneous multislice accelerated sequences may also foster the advent of shorter TRs in liver IVIM studies in the future. 37,38 Thus, a further evaluation of short TR settings might be needed.
We acknowledge several limitations to our study. First, all images were acquired with free breathing only. Respiratory trigger techniques such as prospective acquisition correction could have been used, which would have reduced the effect of liver motion during the image acquisition. Due to this motion, images that should ideally be at the same slice position might have been acquired at different slice positions. Concerning respiratory-triggered examinations for which the difference between adjacent slices is of particular importance owing to saturation effects, interchanged slices might falsify or at least alter the results. In addition, image quality is reduced and, at worst, motion artifacts can occur. 39 Image acquisition with respiratory trigger systems prolongs the acquisition time and does not fix TR, which were the reasons for choosing free breathing in this study.
Additionally, D * was determined using a biexponential IVIM fit. However, this procedure is based on an assumption of limited applicability, as the data showed triexponential behavior at low b values, as shown in previous studies. 5,6,8 We believe that the triexponential behavior originates from a distribution of varying flow velocities of blood and rather represents a data representation than a true biophysical model. This distribution may result from different (blood) compartments and different vessel sizes. 5,40,41 Using diffusion encodings that switch on more than one physical gradient would have been more b-value efficient and would have allowed the sampling of shorter TE. We did not aim for short TE, however, in order to keep the perfusion fraction sufficiently large, 4 and thus used the single physical gradient encoding scheme.
Finally, although the sample size of 15 volunteers is larger than the average sample size of many other MR studies involving healthy volunteers, 42 an even larger sample size might have revealed significant dependencies on the slice setting. Nonetheless, with our measurement precision and setup, we deem such dependencies to be not of major relevance. They may become more relevant if further technical advances improve the precision, or if other setups-for example, shorter TR-are used.

CONCLUSION
In conclusion, this study showed that the biexponential IVIM parameters D, f , and D * are comparable among IVIM studies of the liver that use different slice settings. Saturation effects only had a small effect on the obtained IVIM parameters, including the perfusion fraction, although the mean f slightly decreased in the many slice setting. This small size of saturation effects, however, only holds true for TR that are approximately as long as those used in our study and may be different for shorter TR.

FUNDING INFORMATION
Funding from the Deutsche Forschungsgemeinschaft (DFG) is gratefully acknowledged (DFG project 446875476).

ACKNOWLEDGEMENT
Open Access funding enabled and organized by Projekt DEAL.

SUPPORTING INFORMATION
Additional supporting information may be found in the online version of the article at the publisher's website.
FIGURE S1. Estimation of flow effects in the liver. Grey slices indicate even slice numbers. Blood flow is described by the probability W of blood flowing from even to odd (propagation time of TR/2) and from odd to odd slices (propagation time of TR). The saturated blood recovers according to T 1 -relaxation depending on the recovery time (TR/2 for propagation from even slices, TR for propagation from odd slices). FIGURE S2. Estimation of saturation effects for three different T 1 relaxation times of blood (1900, 1600, and 1300 ms). Modeling was performed assuming Gaussian and laminar isotropic pipe flow, respectively. Saturation was simulated for both the MS setting and the excitation of only one slice (corresponding to FS). The dashed horizontal lines indicate the measured mean f in the MS setting (green) and in the FS setting, corresponding to f ∞ (red). The dotted and dashed blue lines indicate the estimated (within TR) from the MTT and the pseudodiffusion coefficient, respectively. FIGURE S3. Estimation of saturation effects for three different TR (3600 ms, 2500 ms and 2000 ms) at T1 = 1300 ms. Modeling was performed assuming Gaussian and laminar isotropic pipe flow, respectively. In addition, saturation was simulated for both the MS setting and the excitation of only one slice (corresponding to FS). The dashed horizontal lines indicate the measured mean f in the MS setting (green) and in the FS setting, corresponding to f ∞ (red), respectively for TR = 3600 ms. The dotted and dashed blue lines indicate the estimated (within TR) from the MTT and the pseudodiffusion coefficient, respectively. FIGURE S4. Comparison of IVIM fit parameters for arithmetic and geometric averaging over diffusion directions. For each volunteer and evaluation approach, the IVIM parameters of the FS and MS setting were averaged. FIGURE S5. Comparison of boxplots for D * between fitting with and without S 0 . In summary, values for D * and the p-value remain similar. The fit with S 0 suffered less outliers.

APPENDIX A
Appendix A provides a brief outline of an estimate of the size of the saturation effect with a simple biophysical model (Figures 4 and S1-S3). For a detailed description, see the Appendix S1.
The slices are enumerated (Figure 4), and the zeroth slice is excited at t = 0. The blood and liver tissue undergo T 1 relaxation according to 1 − exp after an excitation pulse, where t is the time that passes after the excitation pulse.
After excitation, the blood flow into adjacent slices is modeled by a propagator according to a Gaussian distribution and a propagator describing laminar isotropic pipe flow. A probability W(n, t) for blood translating by n slices in time t can be obtained by integrating over the propagator accordingly.
W is then used to describe blood flow from even to odd (with propagation time TR/2) and from odd to odd slices (propagation time TR). By accounting both blood propagation and the T1 relaxation times (Figure 4), one can estimate the size of the saturation effect.