Comparison of methods for intravoxel incoherent motion parameter estimation in the brain from flow‐compensated and non‐flow‐compensated diffusion‐encoded data

Joint analysis of flow‐compensated (FC) and non‐flow‐compensated (NC) diffusion MRI (dMRI) data has been suggested for increased robustness of intravoxel incoherent motion (IVIM) parameter estimation. For this purpose, a set of methods commonly used or previously found useful for IVIM analysis of dMRI data obtained with conventional diffusion encoding were evaluated in healthy human brain.


INTRODUCTION
The notion of intravoxel incoherent motion (IVIM) provides a means to describe the effects of the incoherent motion of water in living tissue, in particular diffusion and blood microcirculation, on the diffusion MRI (dMRI) signal. 1 By fitting the IVIM model to dMRI data at multiple levels of diffusion encoding, parameters related to both diffusion and perfusion can be estimated from a single dMRI data set.The use of IVIM has shown promise for example in oncological and neurological applications. 2,3Estimation of perfusion-related parameters via IVIM analysis of dMRI data can potentially also serve as an alternative to dynamic susceptibility contrast MRI and thus reduce the need for contrast agent injections. 4t was early recognized that IVIM parameter estimation is a challenging task, in particular in low-perfused tissues like the brain. 5,6As an alternative to a direct non-linear least squares (NLLS) fit, it was suggested to exploit that the signal attenuation from blood microcirculation is more rapid than that from diffusion, meaning that the effect of diffusion can be isolated and estimated at higher b-values, and kept fixed in a second step when perfusion-related parameters are estimated in a so called segmented algorithm. 5,7,82][13] More recently, deep learning-based algorithms have been suggested to provide the robustness of some Bayesian algorithms, while being orders of magnitude faster. 146][17] The reasoning behind this is that multidimensional diffusion encoding can disentangle model parameters and thus possibly enable a more robust parameter estimation. 18hlgren et al. suggested an implementation in which bipolar gradients were used for diffusion encoding where the second bipolar gradient can be flipped relative to the first one to either obtain a flow-compensated (FC) or a non-flow-compensated (NC) diffusion encoding. 15Moulin et al. and Simchick et al. did, on the other hand, use numerically optimized gradient waveforms for diffusion encoding which facilitates multiple levels of sensitivity to flow for a given b-value. 16,17Following the methodological differences between joint IVIM analysis of FC and NC data, and conventional IVIM analysis, conclusions drawn in previous studies regarding optimal choice of methods for IVIM parameter estimation does not necessarily transfer directly.
In this study, methods for IVIM parameter estimation were evaluated regarding the impact on estimated parameter values and repeatability for joint analysis of FC and NC dMRI data in healthy human brain.

MR imaging
Eight healthy subjects (age 21-25 y, five males) were scanned using a 3T Philips MR7700, software release 5.9.0, with a 32-channel head coil (Best, the Netherlands).The study was approved by the Swedish ethical review authority (Dnr 2020-00029) and all subjects signed informed consent.Diffusion-weighted images were acquired with a spin echo EPI pulse sequence utilizing bipolar gradients for diffusion encoding similar to Ahlgren et al. (Figure 1). 15The implementation was built on a prototype pulse sequence enabling use of arbitrary gradient waveforms for diffusion encoding originally developed at Lund University. 19wo scans were acquired, each with seven b-values (0, 5, 10, 20, 30, 100, 200 s/mm 2 ), with the second bipolar diffusion encoding gradient either equal to or flipped relative to the first one, resulting in FC or NC diffusion encoding, respectively (Figure 1), similar to Ahlgren et al. 15 Corresponding flow weighting factors (c-values) for the NC scans were 0, 0.37, 0.53, 0.74, 0.91, 1.66 and 2.35 s/mm, respectively, while they were all zero by design for the FC scans.Diffusion encoding was applied in six directions (sides of a cube), which may reduce the influence of background magnetic field gradients compared with use of only three orthogonal directions. 20Other imaging parameters were: TE = 80 ms, TR = 3700 ms, acquired voxel size 2 × 2 × 4 mm 3 reconstructed to in-plane 1.875 × 1.875 mm 2 , 1 mm slice gap, 29 slices, halfscan factor = 0.79, SENSE = 1.9 (anterior posterior direction), phase encoding bandwidth = 27.2Hz/pixel, fat suppression by SPIR and gradient reversal.Dual readout was used for mitigation of Nyquist ghosting. 21An additional b = 0 image with reversed phase encoding direction was acquired to enable correction for susceptibility-induced distortions.The complete set of scans was repeated without moving the subject, thus resulting in two sets of Schematic illustration of the pulse sequence utilized in the study with actual timing displayed.FC and NC denote flow-compensated and non-flow-compensated, respectively.scans with equal measurement setup, to enable assessment of repeatability.Additionally, a scan with six b = 0 images acquired in direct succession was included for assessment of SNR.A T1-weighted 3D image with full coverage of the brain and voxel size of 0.65 × 0.65 × 1.0 mm 3 (sagittal excitation) was acquired for anatomical reference.

Preprocessing of data
The susceptibility distortion field was estimated with TOPUP (FSL v.6.0.4) using the six repeated b = 0 images from the SNR scan and the b = 0 image with reversed phase encoding direction as input. 22,23The average of the six repeated b = 0 images with TOPUP correction applied served as definition of the dMRI space.The TOPUP correction was applied to the remaining dMRI data using APPLYTOPUP (FSL v.6.0.4) and the script eddy_correct (FSL v.6.0.4) was used for correction of motion and eddy currents.The resulting dMRI data were averaged across diffusion encoding directions using the geometric mean.Maps of SNR were generated by calculation of the ratio of the voxelwise mean and SD of the six repeated TOPUP-corrected b = 0 images in the SNR scan.
A segmentation of white matter, cortical gray matter, and deep gray matter was generated with FSL (v.6.0.4) based on the T1-weighted image resampled to 1 mm isotropic resolution. 24,25Voxels outside the brain were excluded with BET and the remaining voxels were segmented into white matter, gray matter and CSF using FAST. 26,27The gray matter voxels were further divided into cortical gray matter and deep gray matter by segmentation of the deep gray matter with FIRST. 28The final segmentation was transferred to dMRI space based on the transformation matrix obtained from an affine registration of the T1-weighted image to the TOPUP-corrected average b = 0 image using FLIRT. 25

IVIM parameter estimation
IVIM parameters were estimated by fitting the IVIM model suggested by Ahlgren et al. to data: where S(b, c) is the signal acquired with diffusion and flow encoding described by b and c, respectively, S 0 is the signal with b = c = 0, D and D b are the diffusion coefficients of water in tissue and blood, respectively, f is the perfusion fraction and v d is the velocity dispersion of blood. 1,15To stabilize the fit, D b was set to 1.75 μm 2 /ms as suggested by Ahlgren et al. 15 The particular expression for signal attenuation due to incoherent flow used in Eq. 1 builds on an assumption of a Gaussian velocity distribution with zero mean and variance 2v 2 d . 15

Non-linear least squares
A direct NLLS fit of Eq. 1 to data was obtained using the trust region reflective algorithm.Start values for the parameters were D = 1 μm 2 /ms, f = 10%, and v d = 1 mm/s, and S 0 was set to the maximum signal value for each voxel.The parameters were limited to the ranges where max S i is the maximum measured signal in a voxel.

Segmented algorithm
The segmented algorithm was implemented such that NC data with b ≥ 100 s/mm 2 were used to estimate D.
Assuming that v 2 d c 2 ≫ 1, Eq. 1 simplifies to which was used in this first step of the algorithm.D and the intercept term A were fixed in a second step where the remaining parameters were estimated using a specifically designed numerical optimization procedure exploiting analytically derived derivatives as described previously. 8he same parameter ranges were used as for the NLLS fit.

Bayesian parameter estimation
A Bayesian algorithm for parameter estimation, building on Markov chain Monte Carlo (MCMC) using Gibbs sampling and the Metropolis-Hastings algorithm, was implemented as described previously. 9The objective of the algorithm is to produce sampled approximations of the marginal posterior parameter distributions from which parameter estimates can be derived, here by calculating the mean value.The sampling process builds on a random walk through the model parameter space where steps are accepted or rejected based on the change in posterior probability, that is, the product of the likelihood function and the prior distribution, for the particular step.An initial 1000 steps (burn in steps) were used to tune the step length to achieve an acceptance rate of approximately 0.5 followed by another 5000 steps with sampling at each accepted step.The number of steps were chosen empirically to ensure convergence of the algorithm by evaluation of multiple runs on the same data.The same parameter ranges as for the other estimation algorithms were implemented as uniform prior distributions over the specified ranges.
In addition to parameter estimation with uniform prior distributions, parameter estimates were also obtained with a spatial prior distribution, which was given by: where  i is the parameter vector of voxel i,  N i are the parameter vectors in the in-plane four-neighborhood N i around voxel i, and || ⋅ || 1 is the L1 norm. 0 was set to D = 1 μm 2 /ms, f = 5%, v d = 2 mm/s or S 0 = max S i to set the differences on a similar scale.Unlike for example, Spinner et al. the denominator was kept fixed to avoid the potential bias that may be introduced if it is updated dynamically. 11

Deep learning-based estimation
The deep learning-based algorithm suggested by Barbieri et al. 14 and further optimized by Kaandorp et al. 29 was modified to enable parameter estimation based on Eq. 1.
Fully connected networks with two hidden layers and a number of neurons in each layer equal to the number of measurements (unique combinations of b-values and c-values) were designed to output each IVIM parameter, which were fed into Eq. 1 to predict the signal.The loss function used for training was defined as the sum of squared difference between the measured and the predicted signal.The output layer of the network utilized a sigmoid function to limit the IVIM parameters to specific bounds (the same as for NLLS).Data were split into a training set (90%) and a validation set (remaining 10%).
Training of the network used the Adam optimizer with a learning rate of 0.03 and dropout of 10%. 30Training continued until no improvement in the loss function, evaluated on the validation set, was seen for 10 consecutive epochs.The implementation was run in a Python 3.9.18code environment with direct calls to pytorch 2.1.0,numpy 1.26.0,nibabel 5.1.0,and tqdm 4.66.

Simulations
Artificial data were generated for separate evaluation of estimation bias and variability.Based on the parameter maps obtained with the Bayesian algorithm using spatial prior distributions, signal values were calculated with Eq. 1 and Rician noise was added at SNR levels (25, 50, … , 275, 300) for all subjects and scans.Parameter maps were then estimated with all the methods described above in the same way as for the measured data.Additionally, equivalent simulations were run on a digital phantom constructed from a 100 × 100 × 50 array, where typical IVIM parameter values were assigned to concentric regions in the 100 × 100 planes given by squares with consecutively decreasing side of 30 pixel, resulting in four regions corresponding to white matter, cortical gray matter, deep gray matter, and CSF (Figure S1).Simulated IVIM parameters were D = [0.85,0.95, 0.90, 1.70] μm 2 /ms, f = [1.5, 2.5, 2.0, 3.5] %, v d = [1.75,1.75, 1.5, 1.5] mm/s, corresponding to the typical results seen in the in vivo data, and S 0 = 1.

Identification of example data
A typical subject was identified by ranking of difference in parameter estimates and repeatability relative to group medians.For each subject, IVIM parameter, estimation method and region of interest (ROI), the estimated parameter value and the difference between repetitions were ranked based on the distance from their respective group median value.The subject with the overall lowest rank was identified as the example subject.Typical voxel data were found in a two-step procedure for each ROI.First, voxels with NLLS parameter estimates within 15% from the group median were identified.Second, among these voxels the one with the median R 2 of the fit to data was selected as example voxel.

Evaluation of parameter repeatability, bias, and variability
Repeatability was calculated as the voxelwise within-subject SD based on the two repeated sets of scans: where x i is the parameter value at measurement i and x = (x 1 + x 2 )∕2.The in vivo within-tissue parameter variability was quantified as the interquartile range within the segmentation of each tissue type.
The bias and variability of parameter estimates obtained from simulations were calculated as the mean and SD of the voxelwise difference between estimated and simulated parameter value, respectively.

Statistical analysis
Difference in estimated IVIM parameter values and their repeatability among estimation methods were tested with the Friedman test using the exact method to calculate the test statistic with the implementation in the Python package permutation-stats. 31When the Friedman test was statistically significant (p < 0.05), pairwise differences between methods were tested with the signed-rank test.

Software
All processing of data, unless otherwise stated, was implemented in Python 3.

Measured data
Acquired data followed the expected behavior with around one or a few percent lower signal at higher b-values for NC data compared with FC data (see example in Figure 2).SNR levels derived from voxelwise fitting residuals were approximately 100-150 (Table 1), that is, average noise fluctuations were on the same level or smaller than the signal difference between FC and NC data that is intended to be observed.The SNR levels derived from repeated measurements were in general higher than those obtained from fitting residuals and were instead around 200 (Table 1).
Estimated IVIM parameter values differed to various extent between estimation methods and brain regions (Table 2).The estimated diffusion coefficients (D) were on average approximately 0.1 μm 2 /ms higher in cortical gray matter relative to white matter and deep gray matter, and estimates obtained with the Bayesian algorithm with uniform prior distributions were on average approximately 0.1 μm 2 /ms lower than estimates obtained with other methods.The estimated perfusion fractions (f ) were typically higher in gray matter than in white matter, although differences were small.Estimates obtained with the Bayesian algorithm with uniform prior distributions were on average 6%-10% higher than the other methods, whereas the NLLS estimates were lower (on average approximately 1%) and the deep learning-based estimates were higher (on average approximately 1%) than the remaining two methods.The estimated blood velocity dispersions (v d ) did not differ between tissue types in a structured pattern, but NLLS estimates were in general at least a factor two lower.The deep learning-based estimates were, on the other hand, a factor two higher and highly similar across tissue types, even though complementary simulations indicated a correct implementation (Supporting text and Figure S4).The same trends in data with regard to difference between estimation methods could be seen in the example parameter maps given in Figure 3 (upper panel) and when visualizing the distribution of parameter data pooled over all subjects (Figure 4, left column).Additionally, the parameter maps obtained with the Bayesian algorithm with spatial prior distributions and the deep learning-based algorithm were less noisy in their appearance (Figure 3, upper panel).
The repeatability of D and f estimated with the Bayesian algorithm with spatial prior distributions and the deep learning-based algorithm were on a similar level or better when compared with the other methods (Table 3).The repeatability of the diffusion coefficient (D) was better in white matter than in gray matter and better when using NLLS, the Bayesian algorithm with spatial prior distributions or the deep learning-based algorithm relative to the other methods.The repeatability of the perfusion fraction (f ) was a factor two better for the Bayesian algorithm with spatial prior distributions and the deep learning-based algorithm relative to NLLS and the segmented algorithm and about a factor 10 better than the Bayesian algorithm with uniform prior distributions.The repeatability of the blood velocity dispersion (v d ) was of similar magnitude for all tissue types and estimation methods except for the deep learning-based method which on average were an order of magnitude lower.The within-tissue variability followed similar trends (Table 4).In particular, variability of v d was very small for the deep learning-based algorithm.

Simulated data
The parameter maps and values obtained from simulated data agreed well with those obtained from measured data for a matched SNR level, both with regard to visual appearance of parameter maps (Figure 3) and the trends in estimated parameter values between estimation methods (Figures 3 and 4).Estimation bias and variability had similar SNR dependencies for all tissue types (Figures 5 and 6).The SNR dependence of the estimation bias of D and f were similar, but of opposite sign where D was negatively biased and f was positively bias (Figure 5).NLLS and the Bayesian algorithm with uniform prior distributions were associated with a larger estimation bias and strong SNR dependence.Estimates of v d were more strongly biased with the segmented algorithm and the deep learning-based algorithm.Estimation variability of D was lower using the Bayesian algorithm with spatial priors and the deep learning-based algorithm (Figure 6).The same holds for estimates of f except at higher SNR levels where the segmented algorithm performed similar to the two superior algorithms.Conversely, for v d the estimation variability was lowest with the Bayesian algorithm with spatial priors, highest with the segmented algorithm and mainly insensitive to the SNR level for the deep learning-based algorithm.Equivalent results regarding bias and variability could be seen with the digital phantom, although with somewhat more exaggerated trend for v d , for example, with zero variability for the deep learning-based algorithm (Figures S2 and S3).

T A B L E 2
Estimated IVIM parameter values in different brain regions and for different estimation methods.

F I G U R E 3
Intravoxel incoherent motion parameter maps estimated with the different methods investigated in this study for an example subject (same as in Figure 2).The upper section shows maps estimated from measured data, whereas the lower section shows maps estimated from artificial data generated from the parameter maps obtained from the Bayesian algorithm with spatial priors with applied noise at a SNR level of 150, similar to that seen in the measured data.D, diffusion coefficient; f , perfusion fraction; v d , blood velocity dispersion.

DISCUSSION
In this study, five methods for IVIM parameter estimation based on joint analysis of FC and NC dMRI data were evaluated for application in healthy brain.Results both from measured data and simulations show that a Bayesian algorithm with spatial prior distributions was superior to the other estimation methods with regard to measurement repeatability and estimation accuracy, although a deep learning-based algorithm showed high resilience to noise but with an associated tendency to estimation bias.The results of the current study also add to the currently scarce literature on IVIM parameter values from joint analysis of FC and NC dMRI data in healthy human brain. 15ive methods for IVIM parameter estimation were selected for comparison in the current study.2][13] NLLS and segmented algorithms are arguably the most commonly used methods for IVIM parameter estimation, but Bayesian algorithms are commonly found superior.Bayesian algorithms are more of a family of methods with a strong dependence on the choice of prior distribution, 9 two algorithms with distinctly different prior distributions were included in the evaluation where the uniform prior distribution represents non-informative priors and the spatial prior distribution is more of an informative prior.A Bayesian algorithm with a hierarchical prior has also been suggested by Orton et al., which can enhance the quality of IVIM parameter maps, but was not included in the current comparison as it is designed for analysis of specific regions of interest with a relatively limited true range of IVIM parameter values and, therefore, may be less suitable for organs or regions with multiple tissue types. 32Deep learning-based IVIM parameter estimation has received much attention in recent years and is under rapid development. 14,29,334][35][36][37] Given the interest on this topic in the research community, it is likely that Difference in estimated parameter values between repeated scans in different brain regions and for different estimation methods.

Note:
The displayed values are group median (25th percentile, 75th percentile) of ROI median voxelwise absolute differences with superscript index (i) indicating statistically significant difference (p < 0.05) between the particular method and method i.The methods are abbreviated as NLLS = nonlinear least squares (i = 1), Segmented = segmented fitting algorithm (i = 2), B: Uniform = Bayesian algorithm with uniform priors (i = 3), B: Spatial = Bayesian algorithm with spatial priors (i = 4), DL = deep learning-based algorithm (i = 5).IVIM parameter are D = diffusion coefficient, f = perfusion fraction, v d = blood velocity dispersion.

T A B L E 4
Interquartile range of parameter values within different brain regions for different estimation methods.

Note:
The displayed values are group median (25th percentile, 75th percentile) of ROI median values with superscript index (i) indicating statistically significant difference (p < 0.05) between the particular method and method i.The methods are abbreviated as NLLS = nonlinear least squares (i = 1), Segmented = segmented fitting algorithm (i = 2), B: Uniform = Bayesian algorithm with uniform priors (i = 3), B: Spatial = Bayesian algorithm with spatial priors (i = 4), DL = deep learning-based algorithm (i = 5).IVIM parameter are D = diffusion coefficient, f = perfusion fraction, v d = blood velocity dispersion.

F I G U R E 5
Parameter estimation bias as a function of SNR for different estimation methods in different ROIs.For each method, subject, ROI and parameter, the mean difference between the estimated value and the true value was obtained and shown in the plot as black connected markers (group mean) with error bars (group SD).The vertical thin black line shows zero bias.The methods are abbreviated as NLLS, nonlinear least squares; Segm., segmented algorithm; B: Uni., Bayesian algorithm with uniform priors; B: Spat., Bayesian algorithm with spatial priors; DL, deep learning-based algorithm.D, diffusion coefficient; f , perfusion fraction; v d , blood velocity dispersion.
deep learning-based methods will improve substantially in the near future. 33here was a general good agreement between the results from measured data and those from simulated data at matching SNR levels, including similar spatial patterns in parameter maps (Figure 3) and corresponding differences in estimated parameter values among estimation methods (Figure 4), thus suggesting valid simulations.Based on this, it is clear from the simulations that the Bayesian algorithm with spatial prior distributions and the deep learning-based algorithm can be even more superior to the other estimation methods at lower SNR.The matched SNR level was based on SNR estimates calculated from the fitting residuals as this more closely resembles the discrepancies in data that the estimation methods have to process, including not only purely random noise but also semi-random effects like residual gross motion or pulsation.The SNR levels calculated from repeated Parameter estimation variability as a function of SNR for different estimation methods in different ROIs.For each method, subject, ROI and parameter, the SD of the differences between the estimated values and the true values was obtained and shown in the plot as black connected markers (group mean) with error bars (group SD).The methods are abbreviated as NLLS, nonlinear least squares; Segm., segmented algorithm; B: Uni., Bayesian algorithm with uniform priors; B: Spat., Bayesian algorithm with spatial priors; DL = deep learning-based algorithm.D, diffusion coefficient; f , perfusion fraction; v d , blood velocity dispersion.measurements were higher than their equivalent based on fitting residuals.This is expected as the fitting residuals also contain systematic deviations from the model.One can, for example, see in Figure 2 that the model offsets from the measured data at b = 0 in some cases.This may be due to a more tortuous capillary system in the particular tissue which makes the assumption of a ballistic microcirculatory flow regime invalid, 15,38 but may also be due to a partial volume effect with CSF.The difference between the two SNR calculations may also be due to the difference in scan time for the underlying data.Factors like signal drift become apparent at longer scan times and may manifest as increased residuals. 39 Bayesian estimation algorithm with spatial prior distributions appears to be preferable for joint IVIM analysis of FC and NC dMRI data given the results in the current study.This is consistent with previous findings for IVIM with conventional diffusion encoding as well as for other branches of MRI-based parameter mapping. 11,40,41Data, however, should be carefully inspected if focal lesions are of interest as local structures have been found to faint for some Bayesian implementations. 42It is also worth noting that the other Bayesian algorithm evaluated in this study, using uniform prior distributions were among the least favorable methods.This emphasizes the fact that a general statement that Bayesian algorithms are preferable for IVIM analysis cannot be made.Instead, the choice of prior distribution is of particular importance and also possibly to some extent the choice of central tendency measure used to obtain parameter estimates. 9A deep learning-based algorithm showed similar bias and variability as the Bayesian estimation algorithm with spatial prior distributions regarding D and f , but tended to collapse to estimation of a constant, high, v d .No of the other evaluated methods produced maps of v d with particularly high quality, but the diffuse maps generated by for example, the Bayesian estimation algorithm with spatial prior distributions did at least give estimates close to the expected values.Computational time is also a factor that one may need to consider in the selection of estimation method.The segmented fitting algorithm is inherently fast as it reduces the dimensionality of the numerical problem and opens up for efficient, specifically tailored implementations. 8This is in contrast to Bayesian algorithms with spatial prior distributions where the computational demands are substantially higher.In view of this, the segmented fitting algorithm can serve as a rapid first quality check of data, in particular for D and f , while a Bayesian algorithm with spatial prior distributions can provide high-quality parameter maps when computational time is less of an issue.Deep learning-based algorithms may serve as a compromise, being faster than Bayesian algorithms but still being robust to noise.The implementation used in this study trains the networks for each subject, but if pretrained networks can be used a substantial increase in computational speed can be obtained.
The model used for joint analysis of FC and NC data in this study (Eq. 1) makes two important assumptions.First, it assumes that the ballistic regime is valid, that is, that the velocity of blood is constant, both in terms of magnitude and direction, throughout the diffusion encoding.An encoding time dependence has been observed in human liver and pancreas as well as mouse brain for encoding times relevant to clinical MRI. 38,43The potential encoding time dependence in human brain remains to be determined, but the results of the current study and Ahlgren et al both display a separation between FC and NC data indicating that the temporal characteristics at least are close to the ballistic regime for the employed encoding time (here 50 ms).Second, the model assumes a Gaussian velocity distribution.Assuming isotropic pipe flow results in a heavy-tailed velocity distribution where a non-negligible part of the signal from the perfusion compartment remains at b-values where such contributions typically are assumed to be zero (∼100 s/mm 2 ). 44Disentangling this potential residual signal component from the signal of the diffusion compartment to confirm or reject this hypothesis is challenging given the small perfusion fraction and typical noise levels.Violation of the two assumptions would both result in an under estimation of the perfusion fraction, for assumption one due to a too low FC signal and for assumption two due to a too high NC signal.This should be taken into account when comparing results with studies using conventional diffusion encoding.
The pulse sequence and range of b-values were chosen to be similar to those used by Ahlgren et al. 15 In particular, the maximum b-value was 200 s/mm 2 , much lower than the typical maximum b-value for conventional IVIM (around 800-1000 s/mm 2 ).This can be seen as a direct consequence of lower efficiency of bipolar gradient pulses.That is, to not lose SNR, the maximum b-value must be reduced from approximately 800 s/mm 2 to to 200 s/mm 2 to keep the TE constant.However, disentangling of model parameters by separately varying diffusion and flow weighting also opens up for use of lower maximum b-values while still being able to separate the signal components related to diffusion and perfusion.With the high demands on SNR for IVIM parameter estimation in the brain, this can be highly valuable as use of high b-values will either introduce an estimation bias due to "kurtosis effects" which are small at b < 1000 s/mm 2 but relevant at high SNR and small f , or necessitate inclusion of such effects in the model, thus increasing the demands of data even further. 45Using a maximum b-value as low as 200 s/mm 2 has the most direct implications on the segmented fitting algorithm as the threshold b-value was set to 100 s/mm 2 instead of the often-used 200 s/mm 2 . 46iven the shape of observed signal versus b curves, this seems to be a reasonable assumption for healthy brain.The smaller range of b-values used in the first step of the segmented algorithm, however, is a likely explanation for the higher variability of D estimated with the segmented algorithm (Table 3).
In addition to providing additional information for IVIM modeling, FC diffusion encoding has been shown to be of value in multiple other applications.It does, for example, mitigate signal voids in dMRI of the left liver lobe originating from the pulsative motion of the heart, 47,48 enable cardiac dMRI, 48,49 and can reduce phase errors in segmented 3D EPI-based dMRI. 50Use of FC diffusion encoding, however, does remove much or all of the "IVIM effect" in measurement setups close to the ballistic regime, making it non-trivial to combine the improved modeling capabilities with variable flow encoding and the increased robustness of flow-compensation.
The study is associated with some limitations.The study cohort is small although the main results are distinct and consistent across subjects.The study also only used a single set of b-and c-values.At the current stage of development, exploratory protocols like the one used in the current study are warranted, but for larger studies and clinical deployment optimized protocols are more suitable for which specific evaluations of estimation methods may be needed. 17,46,51,52Additionally, the simulations based on in vivo parameter maps used one of the evaluated methods as ground truth, which may bias those results to some extent.A fully artificial digital phantom, for which such bias is not applicable, was also used for simulation, showing similar results, but such phantoms typically have simplistic geometry and, in particular, lack realistic within-tissue variability of IVIM parameter values.Development of realistic digital reference objects would be of particular value.Finally, the deep learning-based method was adapted to the ballistic regime without any further optimization.Complementary simulations showed that the observed inability to accurately estimate v d was related to the IVIM parameter ranges typical for brain tissue as evaluated in the current study, while evaluation on broader parameter ranges equivalent to those used by Kaandorp et al.  resulted in reasonable performance (Supporting text and Figure S4). 29,33A possible explanation could be that with only small values of f , the signal dynamics introduced by varying v d are very subtle and harder for the network to learn.These results warrant further studies developing new network architectures and training paradigms suitable for the IVIM parameter ranges typical for brain tissue.

CONCLUSIONS
Choice of estimation method has a strong influence on IVIM parameter maps obtained from joint analysis of FC and NC dMRI data in healthy human brain, where a Bayesian algorithm with spatial prior distributions was found to be preferable although deep learning-based algorithms show promising results warranting further development.

F I G U R E 2
Illustration of ROIs in white matter, cortical gray matter, and deep gray matter (left column), and signal-versus-b curves (middle and right column) for an example subject.The data are either average signal in the ROI (middle column) or typical data from a single voxel (right column).The parameter values displayed in the middle and right columns were obtained from a NLLS fit to the displayed data.FC and NC denote flow-compensated and non-flow-compensated, respectively.
Estimated SNR levels in different tissue types calculated either as the estimated S 0 divided by the SD of the fit residuals (SNR fit ) or as the mean and SD of six b = 0 images acquired in direct succession (SNR rep ) T A B L E 1 11,12As the