Assessment of intravoxel incoherent motion MRI with an artificial capillary network: analysis of biexponential and phase-distribution models.

PURPOSE: To systematically analyze intravoxel incoherent motion (IVIM) MRI in a perfusable capillary phantom closely matching the geometry of capillary beds in vivo and to compare the validity of the biexponential pseudo-diffusion and the recently introduced phase-distribution IVIM model.
METHODS: IVIM-MRI was performed at 12 different flow rates ( 0.2 ⋯ 2.4 m L / min ) in a capillary phantom using 4 different DW-MRI sequences (2 with monopolar and 2 with flow-compensated diffusion-gradient schemes, with up to 16 b values between 0 and 800 s / mm 2 ). Resulting parameters from the assessed IVIM models were compared to results from optical microscopy.
RESULTS: The acquired data were best described by a static and a flowing compartment modeled by the phase-distribution approach. The estimated signal fraction f of the flowing compartment stayed approximately constant over the applied flow rates, with an average of f = 0.451 ± 0.023 in excellent agreement with optical microscopy ( f = 0.454 ± 0.002 ). The estimated average particle flow speeds v = 0.25 ⋯ 2.7 mm / s showed a highly significant linear correlation to the applied flow. The estimated capillary segment length of approximately 189 u m agreed well with optical microscopy measurements. Using the biexponential model, the signal fraction f was substantially underestimated and displayed a strong dependence on the applied flow rate.
CONCLUSION: The constructed phantom facilitated the detailed investigation of IVIM-MRI methods. The results demonstrate that the phase-distribution method is capable of accurately characterizing fluid flow inside a capillary network. Parameters estimated using the biexponential model, specifically the perfusion fraction f , showed a substantial bias because the model assumptions were not met by the underlying flow pattern.

INTRODUCTION

DW‐MRI1, 2 offers the unique possibility to noninvasively assess the microscopic motion of water molecules inside biological tissue. For free thermal molecular diffusion, the measured DW‐MRI signal as a function of the applied diffusion weighting follows an exponential decay. However, DW‐MRI is not only sensitive to thermal molecular diffusion but also to incoherent motion of water molecules in general, coining the term intravoxel incoherent motion (IVIM) MRI.3 A primary source of such motion in vivo is capillary blood flow, which leads to an additional superimposed decay attributed to the signal of flowing blood in perfused capillaries. The exact shape of the signal decay in IVIM‐MRI depends not only on the speed of blood flow but also on the time until a flowing particle changes its movement direction due to capillary curvature or branching. The signal decay due to perfusion after many directional changes during the diffusion experiment, called the pseudo‐diffusion limit, can be approximated by an exponential function. Considering the presence of water molecules in perfused capillaries as well as extravascular space in a single imaging voxel, the measured signal can then be described as a biexponential decay function of the applied diffusion weighting in the pseudo‐diffusion limit.

This biexponential pseudo‐diffusion model enjoys great popularity because it is elegant and seemingly fits measured data from a great variety of organs very well. However, whether the preconditions for the pseudo‐diffusion model are met is often unclear or ignored. This issue was recently addressed by Wetscherek et al4, 5 by introducing a method that does not assume the pseudo‐diffusion limit but is based on phase distributions of the NMR spin ensemble generated by the simulation of particle pathways. With this approach, it was both shown that the pseudo‐diffusion limit was not reached in common IVIM‐MRI experiments in the liver and pancreas and that the new method also allows estimation of the average speed of the blood flow (i.e., the mean value of the statistical distribution of particle speeds) as well as the characteristic duration until a particle changes its movement direction. Furthermore, picturing the microvasculature as a successive network of straight segments, the average length of a single segment can be calculated via . In contrast, the conventional pseudo‐diffusion model combines all these measures into the pseudo‐diffusion coefficient6 .

Unfortunately, the validation of IVIM‐MRI methods in vivo is inherently difficult. Quantities such as the blood flow velocity and the capillary length are extremely challenging, if not impossible, to determine in vivo without the use of highly invasive methods and can hardly be regulated. To obtain a better understanding how the IVIM parameters relate to flow characteristics, various efforts to mimic blood flow have been undertaken using flow phantoms made of plastic tubes,7, 8 sponges,9 or columns packed with microspheres3, 10, 11 to simulate capillary perfusion. However, these approaches deviate severely from the actual geometry of capillary beds in vivo and hence from the theory that IVIM‐MRI is based on.

The purpose of this study was to systematically analyze IVIM‐MRI in a perfusable capillary phantom closely matching the geometry of capillary beds in vivo and to compare the validity of the biexponential pseudo‐diffusion and phase‐distribution IVIM models.

THEORY

Incoherent motion of water molecules in blood caused, for example, by capillary perfusion during the application of diffusion‐sensitizing gradients results in a distribution of phases and, accordingly, in an attenuation of the measured signal intensity. This led to formulation of the IVIM model by Le Bihan et al.,3 which expresses the signal as a function of the diffusion weighting ( value) in terms of 2 compartments with signal fractions and

Both compartments are subject to thermal diffusion with diffusion coefficient , whereas the compartment with fraction experiences additional signal attenuation by further incoherent motion, specifically capillary perfusion, summarized in the factor (which, in general, may depend not only on the b value but also on the specific shape of the gradient profiles).

The factor also depends on the characteristics of the capillary blood flow, in particular on the flow velocity and the time until a particle changes its movement direction with respect to the duration of the diffusion‐sensitizing period.3 Conceptualizing the capillary bed as segments of length that are interconnected, approximately straight, and randomly (isotropically) oriented, a directional change occurs as a particle traverses from 1 segment to the next.

In the limiting case , called the pseudo‐diffusion limit, can be approximated by an exponential function3 with the pseudo‐diffusion coefficient6 . Note that in the pseudo‐diffusion limit, does not depend on the specific temporal profile of the diffusion gradients.

A framework to analyze the signal attenuation caused by an arbitrary number of directional changes during application of the diffusion‐sensitizing gradients was recently presented by Wetscherek et al.4, 5 The method is based on statistical distributions of normalized phases generated via numerical simulation of particle pathways during the application of diffusion‐sensitizing gradients. The distribution of normalized phases is thereby only dependent on the normalized temporal profile , of the diffusion gradient and the number of directional changes . The signal attenuation for arbitrary combinations of , , , and can then be calculated via4, 5

With this method, it is therefore possible to estimate and explicitly. Furthermore, the average capillary length can be estimated via . Analytic formulas for the normalized phase distributions can only be deduced for and (pseudo‐diffusion limit). Following Wetscherek et al.,5 the distributions for were generated by the simulation of particle pathways.

The phase distribution method also allows for the consideration of a particle speed distribution with mean . This can be achieved by introducing the distribution of (dimensionless) relative speeds with mean (i.e., expectation) value . The distribution of normalized phases depends implicitly on the relative speeds because the width (maximum obtainable phase) as well as the average number of directional changes are linearly proportional to the respective relative speed . This can be accommodated for by adjusting the distribution of normalized phases for the relative speeds as

Finally, the phase distribution can be calculated by integrating the speed‐adjusted phase distributions multiplied by the distribution of relative speeds with respect to :

By replacing in Equation 4 by , the speed corresponds to the mean of the respective speed distribution .

It should be noted that for the parameter estimation in the limiting case , a measurement protocol using a single (not flow‐compensated) temporal diffusion‐gradient profile employed at multiple values suffices. In contrast, to facilitate robust parameter estimation using the phase‐distribution method, a combination of flow‐compensated and not flow‐compensated gradient profiles with varying gradient durations should be utilized.

METHODS

Capillary phantom

As previously proposed by Bellan et al.,12 a 3D microchannel network was constructed using melt‐spun sacrificial sugar fibers embedded in a synthetic resin. The original design was customized following the method described in detail by Gaass et al13 to facilitate the application of controlled flow through the capillary phantom at varying flow rates while performing MRI. The constructed phantom consists of an intricate microchannel network connected to 2 bordering macro channels equipped with Luer‐Lock adapters to feed and drain the network.

A total of 9 microscopy images of the phantom were taken using an optical microscope (DM2500, Leica, Wetzlar, Germany) and subsequently analyzed using ImageJ14 to characterize the geometric structure of the capillary network. As shown in Figure 1, the phantom consists of a highly interconnected capillary system strewn with spherical dilations. All visually discernible capillary segments were manually delineated, yielding their length, diameter, and angular orientation within the imaged plane. The diameters of all discernible spherical dilations were measured as well. Furthermore, the ratio of the volume inside the capillaries to the total network volume was estimated assuming an observational error of the individual measurements of for the error propagation.

Figure 1

Analysis of microscopy images. (A) Exemplary microscopy image of the capillary network. (B) Visualization of the angular distribution of the marked capillaries within the imaged plane. The angular distribution was mapped onto a polar rose plot, creating a direction‐intensity histogram. All angles were mirrored on the horizontal axis because a capillary orientation of corresponds to and vice versa. (C–E) Histograms of measured segment lengths, segment diameters, and dilation diameters, yielding averages of , , and , respectively

The capillary network was flooded with a sodium‐chloride solution (213 g NaCl/L, free diffusion coefficient ) matching the magnetic susceptibility of the matrix epoxy resin to facilitate a homogeneous magnetic field inside the phantom during MRI.

MRI setup

MRI was performed using a 3‐Tesla clinical MRI scanner (Magnetom Skyra, Siemens Healthineers, Erlangen, Germany). Before placing the capillary phantom inside a 16‐channel wrist coil, it was connected to a syringe and filled with NaCl solution via Luer‐Lock adapters and extension tubes. A syringe pump (Standard Infuse/Withdraw PHD2000 Syringe Pump, Harvard Apparatus, Cambridge, MA) was then used to generate precise (accuracy , stated by the vendor) flow at varying rates through the capillary network during MRI. Before each measurement, the entire system was visually inspected for air bubbles and vented if necessary.

Imaging

MRI was initiated using a localizer scan with high in‐plane resolution (spoiled GRE sequence, ; flip angle ; reconstructed voxel size ). DW‐MRI was planned sagittally through the center of the phantom to avoid unwanted signal contributions from the macro channels. Twelve different flow rates ranging from to , increasing in increments of , were applied while performing DW‐MRI experiments. For each flow rate, all sequences listed in Table 1 were employed. The respective sequence ID is referring to a monopolar or flow‐compensated gradient profile followed by the leading edge separation time between two opposing gradient lobes in milliseconds. Sequences FC11 and FC18 employed flow‐compensated diffusion gradients following the scheme . Sequences MP28 and MP171 employed monopolar diffusion gradients following the scheme or , respectively (MP171 utilizes a STEAM15 pulse scheme to circumvent excessive signal loss due to T2‐relaxation). The nominal value sampling scheme was = 0; 5; 10; 15; 20; 30; 45; 60; 80; 110; 150; 200; 280; 400; 600; 800 s/mm2; FC11 sampled a reduced value set up to a maximum b value of 200 s/mm2 to achieve shorter diffusion‐gradient durations. The diffusion‐sensitizing gradients were applied in the read (through‐bore) and phase (vertical) direction. To preclude adverse effects on the measured signal decay, diffusion weighting in the slice direction (left–right) was discarded because it aligns with the direction of the macroscopic flow from inlet to outlet. For parameter estimation, note that the corrected values accounting for the effects of the imaging gradients were used. The corrected b values were calculated using gradient simulations of the respective MRI pulse sequence and are provided in Supporting Information Table S1. Further sequence parameters were ; variable (see Table 1); voxel size matrix size ; FOV ; parallel imaging factor 2 (GRAPPA); and 3 signal averages. All flow rates of a series were measured in 1 session with a total duration of approximately 10 hours. A second series was performed on a different date to assess the repeatability of the estimated parameters.

Table 1

Applied DW‐MRI sequences with varying gradient schemes and gradient timings

ID	Gradient Scheme	b values [s/ mm2]	Δ [ms]a	TE [ms]
FC11	flow‐compensated	0‐200	11.2	58
FC18	flow‐compensated	0‐800	17.8	84
MP28	monopolar	0‐800	28.1	60
MP171	monopolar (STEAM)	5‐800	171.2	47

Abbreviations: FC, flow‐compensated; MP, monopolar.

Δ denotes the leading‐edge separation time between and for flow‐compensated and between and for monopolar gradient schemes.

Speed distribution

Accounting for the particle speed distribution inside the constructed phantom is critical to facilitate accurate assessment of the capillary flow characteristics using the phase‐distribution model.4, 5 Based on previously published results16, 17 on the related topic of fluid flow in fracture networks, the speed distribution inside the constructed capillary network was approximated by the following highly asymmetric and L‐shaped speed distribution:where denotes the Heaviside function. Thus, has essentially the shape of scaled such that the 0 occurs at . For a more detailed discourse on the particle speed distribution inside the constructed phantom, please refer to the online Supporting Information (see Supporting Information Figures S1, S2, and S3).

Parameter estimation

All analyses of the acquired MR images were done using MatLab 9.5 (MathWorks, Natick, MA). For each image, the signal in the capillary network was (arithmetically) averaged inside regions of interest based on thresholding of the least diffusion‐weighted image of the respective sequence. A threshold of more than 30 times the SD of the background noise ensured that only voxels with strong signal from the capillary network were included in the regions of interest. The resulting regions of interest consisted of approximately 30 voxels. The measured signal intensities from diffusion weighting along the read and phase directions were (geometrically) averaged.

The model parameters were estimated via nonlinear least‐squares fitting of the respective model to the measured signal attenuation using the MatLab 9.5 function fminunc.18 The fitting of the phase‐distribution model requires the normalized phase distribution (see Theory section) to be generated a priori. Standard errors for the estimated parameter values were approximated as the square root of the diagonal elements of the covariance matrix19, 20; 95% confidence intervals were calculated as . For more detailed information on the fitting procedure of the phase‐distribution model, please refer to the source code provided by Wetscherek et al.5

Note that the pseudo‐diffusion model was fit solely to data measured with sequence MP28, whereas the parameter estimation using the phase‐distribution model comprises all 4 sequences listed in Table 1. Depending on the respective gradient profile, the matrix material of the constructed phantom restricts the Brownian motion of the water molecules. Therefore, fitting separate (thermal) restricted diffusion coefficients for each sequence is required to account for this effect using the phase‐distribution model.

Statistics

To compare the quality of the signal decay model fits to the measured data, the Akaike information criterion21 and the respective weights22 were calculated. To examine the linear proportionality of the estimated average particle speed to the applied flow rate, a linear regression analysis of the estimated speed versus the applied flow rate was performed. To quantify the agreement between the 2 measurement series, A and B, the coefficient of variation23, 24 () was calculated viawith , , and the number of estimate pairs and .

RESULTS

Figure 1A portrays an exemplary microscopy image of the capillary network. The dissolved sugar structures leave a highly interconnected capillary system. Furthermore, the network is strewn with spherical dilations likely originating from larger sugar crystals and residual bubbles of air. In total, capillary segments and spherical dilations were delineated inside the imaged sections.

The angular orientations of the marked capillaries appear to be uniformly distributed (Figure 1B). Figure 1C thorough 1E depict the histograms of the measured capillary segment diameters, segment lengths, and dilation diameters with averages of , , and , respectively. The ratio of the volume inside the capillaries to the total network volume within the imaged sections was estimated at , with the volume inside the spherical dilations .

Figure 2A shows coronal, axial, and sagittal views through the phantom measured with the gradient echo sequence. The fluid‐filled macro‐channels serving as inlet and outlet demonstrate high SNR. The capillary network in between yields lower signal intensity because the matrix material is taking up most of the volume. Exemplary diffusion‐weighted images (Figure 2B) of the capillary phantom using the MP28 sequence at an applied flow rate of show decaying signal with increasing value. The capillary network appears distorted in the phase‐encoding direction (horizontal in presented image orientation), an artifact common to single‐shot EPI sequences.

Figure 2

MR images of the constructed phantom. (A) Coronal, axial, and sagittal views through the capillary phantom using a FLASH34 sequence. (B) Exemplary diffusion‐weighted images of the capillary phantom using the MP28 sequence at an applied flow rate of

Initially, the measured signal decays (measurement series A) were analyzed using a model comprising of a single, flowing compartment:where is based on the generated normalized phase distributions (see Equation 4). However, as illustrated in Figure 3A through 3C, the model fit is lacking, especially at higher flow rates. The residual sum of squares shows a strong dependency on the applied flow rate, and the residuals appear to display a systematic pattern. The introduction of a second, static compartment to the model, sharing the molecular diffusion coefficient with the flowing compartment,shows a profound improvement of the model fit. Furthermore, the remains largely constant over the whole range of the applied flow rates. Comparing the model fits at an exemplary applied flow rate of , the respective Akaike weights are and , indicating that the 2‐compartment model is more appropriate with very high certainty.

Figure 3

Comparison of the model fits using a single, flowing compartment (1‐comp model) and a model comprised of a flowing and an additional static compartment. (A) Exemplary model fit at an applied flow of (only data for sequence MP28 is displayed). (B) Residuals between the model fit and the measurement data (only data for sequence MP28 is displayed). (C) ‐plotted versus the applied flow rate. (D) Measured signal intensities of all sequences and the respective 2‐compartment phase‐distribution model fit at an applied flow of . FC, flow‐compensated; MP, monopolar; RSS, residual sum of squares

Figure 3D displays measured signal intensities and the respective (2‐compartment) phase‐distribution model fit for all 4 sequences at an exemplary applied flow rate of . Note that although the sequences are colored individually, only a single model fit was performed (i.e., fitting a single value of , , and ), incorporating the data of all 4 sequences. In other words, the displayed fitted curves all correspond to the same value of , , and ; however, their slope is different because the diffusion gradient scheme varies between sequences, and therefore the generated phase‐distributions also differ. In general, the phase‐distribution model fits the measured data well, as demonstrated by the fitted curves closely following the measurement points.

A comprehensive illustration of the obtained results using the phase‐distribution model is presented in Figure 4. All parameter estimates and calculated s are listed in Supporting Information Table S2 as well. At low flow rates, the parameter estimates are associated with large uncertainties, especially and . For flow rates of , the confidence intervals are appreciably small and results from both measurement series are in good agreement, yielding s of , , , and for , , , and , respectively. The estimated fraction of the signal attributed to the flowing compartment yields an average of (mean SD of all flow rates). The estimated average particle speed shows a linear proportionality to the applied rate of flow and increases from to , respectively to , over the course of the measurement series. The linear regression of (measurement series A) versus the applied flow rate yielded an coefficient of determination of with a P value of and an estimated y‐axis intercept of . Analogously, measurement series B yielded , , and a y‐axis intercept of .

Figure 4

Estimated parameters and respective 95% confidence intervals versus the applied rate of flow for 2 measurement series A and B. In general, there is a good agreement between the 2 measurement series. (A) remains largely constant closely matching the ratio of the volume inside the capillaries to the total network volume of determined using optical microscopy (dashed line). (B) The estimated average particle speed shows a linear proportionality to the applied rate of flow with narrow confidence intervals. (C) At low flow rates, is associated with large uncertainties but displays an inverse proportionality at rates of . (D) Consequently, the average capillary segment length calculated via remains approximately constant at flow rates of , with estimates close to the segment length of determined using optical microscopy (dashed line). (E) Estimated thermal diffusion coefficients versus applied rate of flow

The characteristic duration behaves inversely proportional to the applied flow rate. Accordingly, the average capillary segment length remains approximately constant at flow rates of with an average estimated length of and (, flow rate ) for measurement series A and B, respectively.

The estimated thermal diffusion coefficients (Figure 4D) exhibit increasing values with decreasing diffusion times Δ, as well as a slight increase over the course of the measurement series.

Figure 5 illustrates the estimated parameters using the biexponential pseudo‐diffusion model compared to the estimates from the phase‐distribution model. For the biexponential model, the particle speed was calculated via , for which a capillary length of (average length determined with phase‐distribution model at rates of flow ) was used. The biexponential estimates for are substantially lower than the estimates obtained using the phase‐distribution model and increase with the applied flow rate. However, the estimates for show a reasonable agreement between the 2 models. The molecular diffusion is estimated consistently higher when using the biexponential model compared to from the phase‐distribution model fit.

Figure 5

Comparison of parameters estimated using the biexponential versus phase‐distribution model. (A) The perfusion fraction is estimated smaller using the biexponential and appears to increase with increase with the applied flow, whereas (B) the particle speed shows a reasonable agreement. (C) Comparison of molecular diffusion coefficients using as the phase‐distribution estimate

A detailed analysis of the residuals and goodness of fits is presented in Figure 6. Note that only the residuals from sequence MP28 are shown here, although the phase‐distribution model fits are based on all 4 sequences. Most prominent at high flow rates, the phase‐distribution model fit results in decidedly smaller absolute values of the residuals. At small flow rates, the biexponential model is preferred by the Akaike information criterion due to the smaller number of fit parameters (Figure 6C). However, as the flow rate increases, the decidedly reduced values (Figure 6B) obtained using the phase‐distribution model outweighs the larger parameter set, and a strong preference of the phase‐distribution model becomes apparent.

Figure 6

Analysis of the residuals and goodness of fits using the biexponential and phase‐distribution model (only data for sequence MP28 is displayed). (A) Fit residuals as line‐plots color‐coded according to the applied flow rate. (B) values of the respective models versus the applied flow rate. (C) Akaike weights at each flow rate. RSS, residual sum of squares

DISCUSSION

In the course of this work, the capability of conventional pseudo‐diffusion (biexponential) IVIM‐MRI and phase‐distribution IVIM‐MRI to characterize capillary flow was investigated and compared. To this end, a phantom facilitating adjustable fluid flow in a capillary bed at in vivo scales was constructed. Based on extensive DW‐MRI measurement series, the potential of phase‐distribution IVIM‐MRI to accurately estimate flow properties, such as the flow speed as well as information about the capillary geometry, was assessed.

Phase‐distribution model

Utilizing the Akaike information criterion, the acquired data was shown to be best described by a 2‐compartment model consisting of a static and a flowing compartment in which the static compartment is hypothesized to be ascribed to liquid inside spherical dilations. The flowing compartment was fit using the phase‐distribution model. In accordance with the expectation that the relative volumes of multiple compartments should be independent of the applied flow rate, the determined signal fraction attributed to the flowing compartment remained approximately constant over all applied flow rates. The average signal fraction of all flow rates from measurement series A is appreciably close to the ratio of the volume inside the capillaries to the total network volume , estimated at using data from optical microscopy. This supports the notion that the static compartment corresponds to liquid inside the spherical dilations. Comparing to in vivo measurements, a perfusion fraction of is relatively large compared to reported values in well‐perfused organs such as the kidneys25 or the liver.5

The estimated thermal diffusion coefficients increase to some extent with the applied rate of flow; however, the majority of estimates is smaller than the determined of of the freely diffusing NaCl solution. Furthermore, the effects of restricted diffusion come into display: the estimated diffusion coefficients show an inverse relation to the leading‐edge separation duration (“diffusion time”) Δ of the respective gradient profile. The slight increase within each measurement series may be caused by the fact that both compartments (i.e., capillaries and spherical dilations) were modeled with a common diffusion coefficient while actually having restriction effects of varying degrees. Furthermore, adjoining capillaries may still evoke turbulences and drifts in the spherical dilations, which will reflect on the estimated diffusion coefficients.

In a closed system, the average particle speed due to flow is expected to be linearly proportional to the applied rate of flow. The estimated average particle speed inside the capillary phantom fully met this expectation. The calculated coefficient of determination between the applied rate of flow and was for both measurement series, indicating that practically all variation in is predictable from the applied flow rate. Furthermore, there is negligible constant bias because the estimated y‐axis intercepts were close to 0. The estimated average flow speeds, ranging from up to , cover reported values in vivo from various studies. Ivanov et al26 measured the red cell velocity in the brain and the temporalis muscle of rats to be and , respectively. Pawlik et al27 reported a median flow velocity of red cells of in the cerebral cortex of cats. The velocity of capillary blood cells in the cutaneous microvasculature was reported to range between to in patients with venous leg ulcers in a study by Stücker at al.28 Due to the associated difficulties, measured values of the blood flow speed in internal organs have never been reported to our knowledge. However, using their proposed phase‐distribution IVIM method, Wetscherek et al4 estimated the speed of flow to be and for the liver and pancreas, respectively.

At very slow flow rates, the characteristic duration until a particle changes its movement direction could not be determined with satisfactory accuracy, resulting in large standard errors. At flow rates greater than , showed an inverse proportionality to the applied flow rate. Consequently, the hypothesis that the capillary segment length, calculated as , is independent of the applied flow rate is confirmed in this domain. Optical microscopy yielded an average capillary length of , which is slightly below the values determined using DW‐MRI ranging from to for flow rates larger than . A possible source of this bias is the missing depth information (in through‐plane direction) in optical microscopy.

Biexponential model

The measured data was analyzed using the biexponential pseudo‐diffusion IVIM model as well; however, the values increased by factors of up to 5 in comparison to the phase‐distribution model, indicating that the pseudo‐diffusion limit was not reached even at the highest flow rate. At slow flow rates, the reduced parameter set of the biexponential model outweighs the reduction in ; yet, as the flow rate increases, the perfusion effects become more dominant and the Akaike information criterion clearly favors the phase‐distribution model. Using the biexponential model, the signal fraction was estimated to be much smaller ( – ) and displayed a strong increase with the applied flow rate compared to the phase‐distribution model estimate. Yet, in translating the pseudo‐diffusion to the particle speed , the estimates from both models display a good agreement. This is especially interesting with regard to the common notion that is most difficult to estimate accurately, suffering heavily in the presence of signal noise and insufficient value sampling.29

IVIM‐MRI phantom studies

To our knowledge, this work presents the first IVIM‐MRI experiments in a perfusable phantom consisting of a 3D microchannel network similar to in vivo capillary beds, and it was the first time that the phase‐distribution model was investigated using phantom measurements. As the most important property, an IVIM‐MRI phantom must facilitate incoherent particle motion within the scale of a single imaging voxel. In this work, this was realized using sacrificial sugar structures as initially proposed by Bellan et al.12 The constructed phantom allowed for the controlled application of fluid flow at adjustable rates, and optical microscopy revealed a capillary structure with comparable dimensions to in vivo capillary beds. However, there are previously published studies based on IVIM measurements in various otherwise constructed phantoms using the pseudo‐diffusion (biexponential) model.

Le Bihan et al3 constructed a phantom consisting of a column packed with porous microspheres with a theoretical perfusion fraction of . However, the analysis of the measured data using the biexponential IVIM model yielded a significantly lower perfusion fraction of only . The authors attribute this discrepancy to “a too‐low water flow (p5).” This is in agreement with the results in the present work: the signal fraction estimated using the biexponential model was shown to suffer from a negative bias and to be dependent on the flow velocity.

In 2012, Cho et al9 presented IVIM measurements in a phantom using encapsulated sponges allowing to apply varying degrees of pressure. Their estimates for displayed a positive correlation to the applied pressure; however, the authors also state that the signal decay using flow‐compensated diffusion gradients remains largely unaffected by the applied pressure. This suggests that the flowing particles were performing a largely linear motion (straight‐flow limit) within the duration of the diffusion experiment. The applicability of the biexponential IVIM model is therefore questionable, which might explain why the estimated signal fraction was increasing with the applied pressure in some of their measurements.

In a recently published study, Lee et al11 performed IVIM measurements in columns filled with gel beads (Sephadex G‐25; GE Healthcare, Piscataway, NJ) of different sizes at varying rates of flow. A flowing water content fraction was estimated for each bead size by analyzing micro‐CT images of the respective column. The signal fraction mostly underestimated flowing water content, especially at low flow rates, and displayed a strong positive correlation to the applied flow rate. Furthermore, the estimated displayed a counterintuitive behavior by initially decreasing as the applied rate flow was increased. In light of the presented work, these occurrences can be explained by a too‐slow flow that thus does not reach the pseudo‐diffusion limit, and 2) by causing little flow‐induced signal decay.

Limitations

In the present work, the speed distribution for the constructed capillary network was approximated by adapting previously published results on the related topic of fluid flow in fracture networks. A comparison of the signal model fits of the experimental data indicated that the determined particle speed distribution in the capillary network is indeed a reasonable approximation of the true speed distribution. However, the model fits using the phase‐distribution model still display residuals with systematic patterns, indicating a remaining discrepancy between the true underlying particle speed distribution and the idealized mathematical model . To approximate the speed distribution for in vivo microperfusion, further considerations will be necessary because flowing blood acts as a non‐Newtonian fluid. The diameters of smaller capillaries can easily be surpassed by the size of human red blood cells ( 30), which are then “squeezed” through the capillaries separated by segments of plasma, creating a plug‐flow effect.31, 32 Naturally, the heterogeneous composition of blood also influences the flow in larger capillaries and microvessels.33

Unfortunately, there is no ground truth for the average particle speed inside the capillary network, and the estimates had to be validated indirectly by performing a regression analysis. As mentioned above, the highly significant linear correlation with the applied flow rate indicates that practically all variation in is predictable from the applied flow rate, and the small y‐axis intercept shows that there is negligible constant bias. This does leave the estimated average particle speed open to a proportional bias, meaning that the true value could be over‐ or underestimated by a certain fixed percentage.

There are also some limitations regarding the constructed capillary phantom. Because most of the volume is taken up by the matrix material, the SNR of a single voxel is quite low. Consequently, a robust voxel‐wise model fitting to generate parameter maps was not feasible. Furthermore, the construction process of the capillary phantom creates a random network of channels, which precludes the manufacturing of identical standardized phantoms. By keeping variables, such as the humidity and the pressure that is applied to compress the sugar fibers constant, it is possible to construct phantoms with similar properties; however, the exact structure of the capillary network will ultimately be unique.

CONCLUSION

In this work, the capability of conventional pseudo‐diffusion (biexponential) IVIM‐MRI and phase‐distribution IVIM‐MRI to characterize capillary flow was investigated and compared in a controllable fashion using a capillary phantom. To our knowledge, it is the first phantom to mimic the key feature of fluid flow inside a highly interconnected network of randomly oriented channels at the scale of in vivo capillary beds. The results show that the advanced phase‐distribution method is capable of accurately characterizing fluid flow inside a capillary network in a reproducible manner, yielding meaningful and intuitive parameter estimates. At the same time, the importance of correct model assumptions emerged because parameters estimated using the biexponential pseudo‐diffusion model showed a significant bias. Specifically, using the biexponential model the perfusion fraction displayed a strong negative bias as well as a dependence on the applied rate of flow.

FIGURE S1 Illustration of fully developed laminar flow profile in a cylindrical capillary (panel A) and the corresponding uniform speed distribution (panel B). The average velocity through a cross section is given by the half of the maximum particle velocity:

FIGURE S2 Histogram of measured capillary lengths (A) and apertures (B) with fitted log‐normal distributions

FIGURE S3 Comparison of the PDFs (solid lines) and CDFs (dotted lines) of the average particle speeds within capillaries in panel A and the overall particle speeds (including also the laminar flow profiles) in panel B for varying shape parameters . The PDFs were normalized to yield an equal expectation value and , respectively

TABLE S1 Nominal and corrected b‐values for each sequence used for the model fit. The corrected b‐values account for the imaging gradients and were obtained through the simulation of the respective MRI pulse sequence

TABLE S2 Parameter estimates for measurement series A and B and resulting coefficients of variation (CV) between the measurement series

Click here for additional data file.