Echo time dependence of biexponential and triexponential intravoxel incoherent motion parameters in the liver

Intravoxel incoherent motion (IVIM) studies are performed with different acquisition protocols. Comparing them requires knowledge of echo time (TE) dependencies. The TE‐dependence of the biexponential perfusion fraction f is well‐documented, unlike that of its triexponential counterparts f1 and f2 and the biexponential and triexponential pseudodiffusion coefficients D*, D1∗ , and D2∗ . The purpose was to investigate the TE‐dependence of these parameters and to check whether the triexponential pseudodiffusion compartments are associated with arterial and venous blood.


| INTRODUCTION
The intravoxel incoherent motion (IVIM) model described by Le Bihan et al. in 1986 treats the measured diffusionweighted signal as a two-compartment model consisting of a diffusion and a perfusion compartment. 1 IVIM imaging has proved useful for detecting and staging various pathologies in many parts of the human body, such as the brain and abdominal organs. 2 In the biexponential IVIM model, the signal decay as a function of the diffusion weighting b is described by a biexponential function, representing the summed signal of perfusion and diffusion compartments: S 0 denotes the unweighted signal strength, f the perfusion signal fraction, D the tissue diffusion coefficient, and D * the pseudodiffusion coefficient.
However, it has been shown that the biexponential signal representation is not well-suited for describing the signal decay at small b-values in some organs, such as the liver [3][4][5][6] and the kidney. [7][8][9][10][11] The observation of a very steep signal decay at very small diffusion weightings made it necessary to reformulate the biexponential signal representation into a triexponential one, the additional "compartment" describing a "faster" perfusion: Here, the term describing the perfusion has been replaced by the sum of two perfusion terms, described by their respective signal fractions f 1 and f 2 and their pseudodiffusion coefficients D * 1 and D * 2 . No consensus has been reached regarding the tissue types and physiological functions to which these two compartments correspond. 5,7 We follow the notion outlined by Riexinger et al, 5 and consider the two perfusion terms in Equation (2) as a data representation 12 rather than representing a biophysical model.
The determination of the IVIM parameters f and, in particular, D * was found to be challenging. [13][14][15] The high uncertainty in the pseudodiffusion coefficient results in a reduced discriminatory value, 16,17 eg, in the differentiation of subtypes of renal neoplasms, 18 or between hepatocellular carcinoma (HCC) and focal nodular hyperplasia (FNH). 19 Despite these difficulties, significant differences in the pseudodiffusion coefficient between pathologies were reported. 20,21 For example, benign, intermediate, and malignant solid softtissue tumors were reported to have significantly different pseudodiffusion coefficients. 22 Further, D * is significantly decreased in fibrous meningioma compared to other meningioma. 23 Recently, the triexponential IVIM approach has also been used for disease characterization. 8,9 IVIM studies of the liver have been performed with a variety of echo times (TEs). For example, the following values have been used: 55 ms, 4 72 ms, 24 100 ms, 5 120 ms. 25 Regarding this range of values, the question arises whether there is a dependency of the IVIM parameters on the choice of TE. If this was the case, the discriminatory value of IVIM parameters would have to be, if not even questioned, at least reassessed in consideration of the used TE. It is already wellknown that the biexponential perfusion fraction f shows a strong dependency on TE in the liver 26 and pancreas, 27 which are organs of relatively short T 2 time. However, besides a conference report, 28 we are not aware of a study that evaluated whether such a TE dependency exists for the other biand triexponential IVIM parameters by acquiring data with different TEs but otherwise fixed settings.
In light of the above described increasing body of evidence that examinations of pseudodiffusion can reveal important information, 22,29 and anticipating that this also holds true for triexponential IVIM parameters, our study was aimed at investigating the dependency of all biexponential and triexponential IVIM parameters on TE, intentionally without using a T 2 correction in the IVIM equation, 26,27 in order to detect the influence on the mere IVIM parameters as they are used in other studies.
A secondary aim is to examine whether the two triexponential perfusion compartments could represent arterial and venous blood by considering their different T 2 times, 30,31 which would alter their signal contribution depending on TE. Given the above-mentioned difficulties in fitting all IVIM parameters, we focused on the liver because it can be investigated well given its size and its large perfusion fraction. 15,27 (1)

| Data acquisition
All data were acquired with a 3T scanner (MAGNETOM Prisma, Siemens Healthineers, Erlangen, Germany) using an in-house developed single refocused spin-echo echo-planar imaging diffusion sequence. 25,32 An 18-channel body coil and the built-in spine coil were used. Guided by an optimization for triexponential IVIM in the liver, 14 diffusion-weighted images were recorded at the following 24 b-values: 0.2 s/mm 2 (two repetitions), 0.3 s/mm 2 , 0.4 s/mm 2 , 0.6 s/mm 2 , 1 s/mm 2 , 1.5 s/mm 2 , 2 s/mm 2 , 3.5 s/mm 2 , 5 s/mm 2 , 6 s/mm 2 (two repetitions), 10 s/mm 2 , 25 s/mm 2 , 35 s/mm 2 , 45 s/mm 2 , 60 s/mm 2 , 70 s/mm 2 , 80 s/mm 2 , 200 s/mm 2 (two repetitions), and 800 s/mm 2 (three repetitions). A relatively large number of small b-values was used in order to correctly sample the initial steep signal decay (c.f. Supporting Information Figure  S1 of Ref. [14]). These b-values were calculated numerically based on the gradient timetable, taking into account imaging gradient pulses. 33 For each b-value, six different diffusion encoding directions were used; (1,1,0), (−1,1,0), (0,1,1), (0, −1,1), (1,0,1), and (1,0, −1), specified in the scanner coordinate system. The numerically calculated b-values differed slightly between the six diffusion directions; therefore, the exact b-value of each diffusion direction was used in the evaluation. This means that six data points with slightly different b-values were used instead of a single data point at the respective nominal b-value. These data were acquired for TE = 45 ms, 60 ms, 75 ms, and 90 ms. The time Δ between the onsets of the two gradient pulses and the diffusion gradient pulse duration δ were kept constant at 21.64 ms and 12.24 ms, respectively, for all TEs.
For the phantom measurement, a vendor-provided bottle phantom filled with saline solution (per 1000 g H 2 O: 3.75 g NiSO 4 × 6 H 2 O, 5 g NaCl) was used. The same protocol was used for the phantom and the volunteer measurements. The temperature was measured with a BOSCH PTD 1 infrared thermometer (Robert Bosch Power Tools GmbH, Stuttgart, Germany). The aim of the phantom measurement was to perform a quality check of the in-house developed sequence.
Fifteen healthy volunteers (age range: 19-58 y, mean age: 24.7 y, median age: 23 y, male/female: 7/8, no known history of liver diseases) were scanned in supine position. Six transversal slices with 4-mm thickness and 4-mm spacing between slices were placed in the liver. Triggering was done in expiration using an external respiratory trigger. Other sequence parameters were field of view (FOV) 400 mm × 400 mm, matrix size 100 × 100 (interpolated to 200 × 200), GRAPPA (acceleration factor 2, 24 reference lines), phase direction anterior-posterior, phase partial Fourier factor 6/8, fat saturation mode: gradient reversal and spectral attenuated inversion recovery (SPAIR), acquisition bandwidth 2272 Hz/ Px, echo spacing 0.5 ms, repetition time (TR) = 1 respiratory cycle. The vendor-provided pre-scan normalize option was used to compensate for surface coil flare. The total acquisition time depended on the volunteers' respiratory frequency and was about 60 min for all TEs and b-values. The study was approved by the Institutional Ethics Committee, and written informed consent was obtained from all volunteers.

| Determination of the IVIM parameters
Data analysis was performed using in-house software developed with MATLAB R2017b (The MathWorks Inc., Natick, MA, USA). The signals of each TE were evaluated separately. For each slice of each TE, a single region of interest (ROI) excluding large vessels was defined on a b = 35 s/mm 2 image and afterward verified, and if necessary, corrected on a b = 0.2 s/mm 2 image. This procedure was chosen because we found it easier to identify and exclude large vessels when blood appeared dark. The left liver lobe was omitted to avoid bias due to the cardiac pulsation artifact, which is more pronounced in the left liver lobe than in the right liver lobe. [34][35][36] Then, the ROI was copied to the images acquired with the other diffusion directions and b-values. The median signal in the ROI was calculated for each acquired image. The median instead of the mean was used to render the evaluation more stable and robust against outliers. These thus calculated signals were normalized to the signal of the lowest b-value. Because the lowest b-value was measured twice and with 6 diffusion encoding directions each, its signal was calculated as the mean signal of all 12 measurements. The normalized data of all slices of each volunteer were used together to fit the IVIM parameters to them, treating the signals of all acquired images (ie, all used diffusion directions for all diffusion weightings) as equally weighted data points. Because S 0 was also fitted (see below), the main purpose of the normalization was to allow a straightforward comparison of fitted curves.
For the phantom measurement, the diffusion coefficient was determined with a monoexponential fit function. The provided error of the diffusion coefficient was the standard error, which was given by the fit algorithm.
For the volunteer measurement, the fit was performed in two steps. A nonlinear least squares fit using the trustregion algorithm was used for each step. The objective function was the sum of squared errors. In the first step, the equation S (b) = S 0 ⋅ exp ( − bD) was fitted to the signals at b = 200 s/mm 2 and b = 800 s/mm 2 to obtain the diffusion coefficient D. The intersection of this monoexponential curve with the signal axis was used to obtain 1 −f , where f is an estimate for the perfusion fraction, which was used to define the bounds. For both the biexponential and triexponential fit, the fit bounds for f, f 1 , and f 2 were set to 0 as the minimum and to 1.5 times the estimated perfusion fraction f as the maximum. The bounds for D * , D * 1 , and D * 2 were set to (0.01, 20) mm 2 /s, (0.005, 0.3) mm 2 /s, and (0.15, 20) mm 2 /s, respectively. The intervals were chosen quite broadly to minimize the influence on the fit result and to avoid the parameters being equal to a bound. In the second step, these bounds were used to fit S 0 and all IVIM parameters except for D, which was kept constant. As it is not possible to measure with an exact diffusion weighting of zero, 33 it was necessary to also fit the initial signal strength S 0 instead of setting S 0 = S (b = 0). To obtain the most accurate set of fit parameters, the biexponential and triexponential fit were each performed 100 times with different starting points, which were randomly selected in the range between the bounds, and the respective set of fit parameters with the lowest sum-of-squares error was saved. No spatial regularization or regularization along the b-value dimension was used. 8 For the biexponential fit, data at b-values between 0.3 s/mm 2 and 6 s/mm 2 were not used in order to maintain consistency with previous studies. 37,38 Each data point was weighted equally, which roughly corresponds to an arithmetic averaging. In an additional evaluation, the signals were averaged geometrically over different diffusion directions (resulting in one signal per b-value) and the same fit was performed. The reliability of the fit procedure was evaluated with a bootstrap approach. 39 For a ROI with n voxels, we sampled the new set of voxels by drawing with replacement n times from the original set of voxels. This means that the new set had also n voxels, with some of the original voxels missing and some of the original voxels being present multiple times. Afterward, the new set of voxels was evaluated analogously to the original one. This resampling (ie, drawing with replacement) was done 1000 times for each ROI, therefore resulting in 1000 sets of IVIM parameters per volunteer and echo time. The ratio of SD, eg, σ f , of the 1000 IVIM parameters and the original parameter value, eg, f, (without resampling) was calculated to assess the reliability.

| Statistics
Statistical analysis was performed with the Kruskal-Wallis test. This test was used instead of one-way analysis of variance (ANOVA) because a preliminary analysis with the Shapiro-Wilk test showed that a normal distribution could not be assumed for any of the measured parameters. For the post-hoc test, the Dunn-Sidak approach was used to correct the p-value for multiple comparisons. For correlation analysis, we calculated the Pearson correlation coefficient. The significance level for all statistical tests was set to 0.05.

| T 2 analysis
The two observed pseudodiffusion compartments can possibly be explained by the venous and arterial blood pool. To evaluate this hypothesis, the expected ratio of their signal fractions f A /f V was estimated by and was compared to the ratio of f 2 and f 1 . The used values for the relaxation times of arterial and venous blood were T 2,arterial = 44 ms and T 2,venous = 107 ms. 31 Figure 1 shows the results of the phantom measurement at TE = 45 ms. The measured diffusion coefficient was D = (2.0620 ± 0.0008) · 10 −3 mm 2 /s. The values for D at 60 ms, 75 ms, and 90 ms were (2.0633 ± 0.0008) · 10 −3 mm 2 /s, (2.0655 ± 0.0008) · 10 −3 mm 2 /s, and (2.0698 ± 0.0007) · 10 −3 mm 2 /s, respectively. For all TEs, the measured signals were very close to the expected monoexponential signal decay. The temperature of the phantom was 21.4°C and 22.2°C before and after the measurement, respectively. The expected water diffusion coefficient for these temperatures is between 2.109· 10 −3 mm 2 /s and 2.153 · 10 −3 mm 2 /s. 40 Figure 2 shows representative images acquired at different b-values and TEs. The ROIs for this volunteer were plotted sparing major vessels and the left liver lobe. As the ROI shape was adapted for each TE, a slightly different shape of the ROI at TE = 45 ms and TE = 90 ms can be perceived.  Figure 3 shows the representative measured signals, bi-and triexponential fit curves, and a comparison of the absolute signal curves at the four TEs for one volunteer. The markers that were used only in the triexponential fit are plotted with lower brightness. Both fit curves match the measured data well; only the biexponential fit curve did naturally not match the excluded data points. The signal decreased at longer TE ( Figure 3B), but the percentage of the initial signal decay at low b-values became larger, indicating the increased perfusion fractions at longer TE.   Figure 4 shows the averaged histograms of the parameter distributions. As mentioned above, for both D * 2 and D * , one of the 60 (15 volunteers, four TEs) calculated relative errors was quite large. The histograms of the distributions leading to these high values are shown as supplemental material (Supporting Information Figure S1, which is available online). Figure 5 shows the resulting IVIM parameters of the volunteer measurements. For each parameter, box plots at the four TEs are shown. Additionally, Table 1 lists the median values (same underlying data as in Figure 5) and results of the statistical analysis with the Kruskal-Wallis test.

| IVIM parameters
A significant dependency of D, D * , D * 1 , and D * 2 was not observed. Outliers were more prevalent for the pseudodiffusion coefficients than for the other parameters; in particular for the triexponential pseudodiffusion coefficients. The "slow" triexponential pseudodiffusion coefficient D * 1 was smaller than the biexponential pseudodiffusion coefficient D * by approximately a factor of two, whereas D * 2 was larger than D * by approximately a factor of 10.
Both the biexponential and triexponential perfusion fractions were significantly dependent on the TE and increased similarly and approximately by a factor of 2 at the largest TE compared to the smallest TE. Table 2 shows the results of the post-hoc test, which was performed for parameters that showed TE dependence. The results for the geometric averaging were almost identical and are presented in the supplemental material (Supporting Information Figure S2).

| T 2 analysis
The calculated values for f A ∕f V and the comparison to the ratio of f 2 and f 1 are shown in Table 3. The values for f A ∕f V decrease with increasing TE, the values for f 2 /f 1 are nearly constant.

| DISCUSSION
In this study, we examined the TE dependence of the biexponential and triexponential IVIM model. The triexponential perfusion fractions f 1 and f 2 and the biexponential perfusion fraction f depended significantly on the TE, whereas the diffusion coefficient D, the triexponential pseudodiffusion coefficients D * 1 and D * 2 , and the biexponential pseudodiffusion coefficient D * did not.
A meta-analysis found that the diffusion coefficient D was (10.9 ± 1.7) · 10 −4 mm 2 /s if averaged over 27 different healthy liver studies. 29 In our study, it lay between 9.4 · 10 −4 mm 2 /s and 9.7 · 10 −4 mm 2 /s, which is within the margin of error of this literature value.
In the same review, 29 the perfusion fraction f was given as 0.231 ± 0.085 and the pseudodiffusion coefficient D * as (7.0 ± 3.1) · 10 −2 mm 2 /s, which is also in line with our results. Table 4 shows an overview of published triexponential IVIM parameters and of those determined in the present study, a graphical representation is shown in Figure 6. The triexponential perfusion fractions f 1 and f 2 were approximately of equal size in all studies with similar TE. The finding that both values increase with TE is confirmed when the previous reports are considered. One exception to this rule is the study by Riexinger et al 5 (cf. TE = 100 ms in Figure 6, see below for a further discussion on the values reported in that study).
The value for D * 2 was generally much larger than that for D * 1 . Even the minimal reported ratio D * 2 ∕D * 1 was > 6 (by Wurnig et al 7 ). In all other studies, D * 2 was larger than D * 1 by at least an order of magnitude, which itself was again at least an order of magnitude larger than D. If the average over all reported D, D * 1 , and D * 2 values is taken, then D * 1 is approximately 25 times larger than D; and D * 2 is approximately 25 times larger than D * 1 . Our finding that neither D nor the triexponential pseudodiffusion coefficients depended significantly on TE is supported when taking into account previous reports; a dependency on TE is not observable in Table 4. The values for D * 1 and D * 2 seem to be conspicuously high at TE = 60 ms, but the values of the quartiles indicate and the statistical analysis shows that this is no significant effect. Our TE-averaged D * 2 value was larger than those of most other studies, which is presumably caused by the lower b-values used in our study. Interestingly, the D * 2 values reported in a previous study 5 were even larger. In comparison to our study, one main difference can be identified: In the present study, we used a higher maximal b-value of 800 s/mm² (vs. 500 s/ mm²). Even though it is often assumed that the kurtosis effect can be neglected for b-values up to 800 s/mm 2 in abdominal F I G U R E 4 Histograms of the parameter distributions obtained with the bootstrap approach. The value on the horizontal axis denotes the ratio between the result of the bootstrap fit and the result of the original fit; that is, a value of 1 means that they are identical.  organs, 7 it might nonetheless have some influence. 41 If this was the case, the signal at high b-values would be larger than expected, leading to a reduced measured diffusion coefficient compared with that at a maximal b-value of 500 s/mm², and may explain the smaller diffusion coefficient in the present study (D ≈ 0.96 · 10 −3 mm²/s) than in the earlier report (D ≈ 1.22 · 10 −3 mm²/s). As a consequence of this smaller D value, the extrapolation of the diffusion compartment to the point b = 0 will hit the signal axis at a lower point, leading to an increased estimate of the perfusion fractions, which could also have an influence on the pseudodiffusion coefficients. This explanation is supported by the finding that the results of another study by the same researcher 14 and a highest b-value of 800 s/mm 2 are more in line with our results. Elucidating the dominant mechanism at work here is beyond the scope of the present study and highlights the challenges associated with obtaining quantitative triexponential IVIM parameters. These challenges are highlighted to an even larger extent if one considers that even the values reported for D, which can be estimated most precisely among the IVIM parameters, 14 vary by 30% in the different studies. The exact same handling of all evaluation steps (ROI placement strategy, fitting process, etc.) appears mandatory if true comparability among studies is sought. Moreover, IVIM parameters may differ significantly between MR scanners from different vendors, 42 which might also be a source of the observed parameter deviations. The evaluation of the fit reliability with the bootstrap approach showed nevertheless that within our study, the averaged relative deviation is smaller than 6% (one outlier left out for D * and D * 2 ), indicating that almost all fit results can be considered trustworthy and stable.
The observation that the perfusion fractions f 1 and f 2 increased with the TE is not unexpected; a similar behavior was observed for the perfusion fraction f in the biexponential IVIM model. 26,27 This can be explained by the different T 2 relaxation times. In the case of the biexponential model, the line of argument is as follows: Blood and liver tissue, which represent the compartments, have different T 2 times, which leads to a TE-dependent value of f. With rising TE, the shorter T 2 time of liver tissue makes the signal fraction of the tissue compartment decrease, and an increase of f can be observed. For the triexponential model, the parameters f 1 and f 2 might be represented by different blood compartments and, therefore, may exhibit different T 2 times. Wurnig et al 7 hypothesized that the compartments might represent arterial and venous blood. Essentially, this would render f 1 and f 2 dependent on TE in a non-identical manner. At 3T, the T 2 time of venous blood is approximately 44 ms; that of arterial blood is approximately 107 ms. 5,31 As T 2,liver was reported to be 34 ms, 43 the perfusion fractions would also increase at larger TEs, but at a very different rate. We observed, however, that the rate of increase was almost equal for f 1 and f 2 , which leads to the conclusion that both compartments have similar T 2 times. The qualitative evaluation strengthens this finding: The ratio of f 1 and f 2 is nearly constant, whereas the expected ratio decreases by a factor of approximately 2 from lowest to highest TE measured. Note that this evaluation only focusses on the general trend of the expected ratio rather than its exact value. For a more exact quantification, the volume fractions of arterial and venous blood would have to be taken into account. However, these are challenging to determine and may depend on the digestive state. 44 The general trend allows the presumption that the fast and slow perfusion compartments contain both arterial and venous blood. This agrees with the observation that a field strength dependency of the bi-and triexponential IVIM parameters does not exist, 5 although venous and arterial blood exhibited different field strength dependencies in terms of their T 2 times. It is important to consider the behavior of f 1 and f 2 when the triexponential parameters are used for disease characterization. Absolute values of f 1 and f 2 can only be a reliable source of diagnostic information when the value of TE is taken into account. If different TEs are used, a bias in f ensues if one compares the values directly. This problem could be overcome by defining a standard TE or some form of TEdependent standard values for f 1 and f 2 , but this appears quite complicated for clinical work. Another approach might be a TE correction strategy as proposed for the biexponential f. 26 The development, application, or evaluation of such a correction strategy, however, was beyond the scope of this work. In contrast, the finding that D * was not dependent on TE indicates that it is possible to compare this parameter among studies performed with different TEs, which is a notable advantage. The same holds true for the triexponential pseudodiffusion coefficients D * 1 and D * 2 . A correction strategy for these parameters does therefore not seem necessary. However, this must be taken with a grain of salt: The diffusion time in our study was kept constant -an important strategy in order to only detect differences due to changes in TE. For example, this avoided any influence of the diffusion time on the diffusion coefficient D 45,46 and on the pseudodiffusion coefficient whose flow correlation time has been reported to be in the order of 100 ms. 25 However, many previous studies coupled TE and the duration of the diffusion encoding. A dependence of D * 1 on the time separation Δ of the gradient lobes was, for example, reported by Fournet et al. 47 Although these measurements were performed in the rat brain, the finding might generalize to further organs and to pathological tissue. For that reason, the use of different TE might still pose some problems concerning the comparability of pseudodiffusion coefficients among studies. Moreover, the large uncertainty associated with the pseudodiffusion parameters should be kept in mind. Potentially, technical improvements will make it possible to measure them with higher precision in the future.
We must also acknowledge several limitations of this study. First, the in-house developed sequence used did not compensate for eddy currents that can induce image distortions. 48 The volunteers who were scanned stated that they had no history of liver disease. However, the images were only checked for obvious pathologies and not by a radiologist, and pre-existing conditions could have falsified the outcome. Additionally, the volunteers were not obliged to follow a prescribed fasting protocol prior to the measurement; different digestive states might possibly have altered the results. Regarding the ROI placement strategy, the use of smaller ROIs might have reduced the risk to accidentally include vessels, but at the expense of statistical power that was obtained with the larger ROIs used. We tried to minimize a possible influence of undesired very bright or dark voxels by the use of the median for ROI signal evaluation. The use of large ROIs may also lead to averaging regions with different diffusion behavior, in contrast to voxel-wise evaluation, where different regions could possibly have been detected. As there is no general consent about the best fit procedure for IVIM, 6 a segmented fit was used and considered the most exact, while other fitting algorithms could have yielded different results. Barbieri et al 49 reported significant differences in the fit parameters for all upper abdominal organs when they compared biexponential fitting algorithms. A further possible limitation (see also above) is that we kept the diffusion encoding pulses identical at all TE, while many vendor-provided sequences adapt their duration to that of TE. This may have limited the comparability to other studies, but the advantage is that the TE dependency can be assessed unequivocally. The dependency of IVIM parameters on the timing of the diffusion encoding gradients remains to be investigated in future studies. Our sample size was limited (15 volunteers), and a larger sample size might have revealed statistically significant TE dependencies for the pseudodiffusion compartments. Nonetheless, it was larger than the reported median of six for MRI volunteer studies. 50 In conclusion, the measured IVIM parameters f, f 1 , and f 2 showed significant dependence on TE, whereas the other IVIM parameters (D, D * , D * 1 , D * 2 ) did not. This indicates that the other IVIM parameters can be compared between studies with different TEs without correction. No indication was found that the two triexponential perfusion compartments represent venous and arterial blood.