Analytical and Numerical Connections between Fractional Fickian and Intravoxel Incoherent Motion Models of Diffusion MRI

Abstract

Impaired tissue perfusion underlies many chronic disease states and aging. Diffusion-weighted imaging (DWI) is a noninvasive MRI technique that has been widely used to characterize tissue perfusion. Parametric models based on DWI measurements can characterize microvascular perfusion modulated by functional and microstructural alterations in the skeletal muscle. The intravoxel incoherent motion (IVIM) model uses a biexponential form to quantify the incoherent motion of water molecules in the microvasculature at low b-values of DWI measurements. The fractional Fickian diffusion (FFD) model is a parsimonious representation of anomalous superdiffusion that uses the stretched exponential form and can be used to quantify the microvascular volume of skeletal muscle. Both models are established measures of perfusion based on DWI, and the prognostic value of model parameters for identifying pathophysiological processes has been studied. Although the mathematical properties of individual models have been previously reported, quantitative connections between IVIM and FFD models have not been examined. This work provides a mathematical framework for obtaining a direct, one-way transformation of the parameters of the stretched exponential model to those of the biexponential model. Numerical simulations are implemented, and the results corroborate analytical results. Additionally, analysis of in vivo DWI measurements in skeletal muscle using both biexponential and stretched exponential models is shown and compared with analytical and numerical models. These results demonstrate the difficulty of model selection based on goodness of fit to experimental data. This analysis provides a framework for better interpreting and harmonizing perfusion parameters from experimental results using these two different models.

1. Introduction

Microvascular perfusion in skeletal muscle is essential for nutrient supply and cellular metabolism [1]. Compromised skeletal muscle perfusion underlies many disease conditions [2,3] and the reduction and loss of mobility with aging [4]. Impaired skeletal muscle perfusion can be related to physical properties such as loss of microvascular function, reduced blood flow, and the rarefaction of microvessels. Diffusion-weighted magnetic resonance imaging (DW-MRI) is a noninvasive measurement that is sensitive to molecular motions in water occurring at different rates within microvasculature and tissue parenchyma, respectively [5]. Thus, the consideration of appropriate signal models and DW-MRI acquisition schemes for quantifying physical properties of these tissue compartments is of particular interest for assessing perfusion in skeletal muscle.

The intravoxel incoherent motion (IVIM) [6] and fractional Fickian diffusion (FFD) models [7] are well-established mathematical representations used to characterize DW-MRI measurements. Compared with the classical single compartmental model which characterizes the random Brownian motion of all tissue water molecules as a monoexponential decay of MRI signal, the IVIM model uses a biexponential representation to differentiate the incoherent motion of intra-vascular blood perfusion and extra-vascular water proton diffusion. The FFD model is a parsimonious representation of anomalous superdiffusion that models the diffusion decay signal as a stretched exponential. In this model, a single diffusion compartment is represented for an exponent value of 1 and superdiffusion is represented when the exponent decreases to values less than 1. This fractional-order model is convenient as it is more robust, with a limited number of diffusion weightings, and more reliable in the presence of noise [8]. However, the relationship between the FFD model parameters and physical tissue parameters is not as straightforward as that of the IVIM model. Comparisons of the biexponential and stretched exponential models have been made extensively in terms of their applications for the diagnosis and grading of disease stages [9,10,11]. The stretched exponential is reported to outperform the biexponential model with a lower estimate variance and higher reproducibility [8,9,10,12]. However, the underlying connection between the two models has not been previously established, making it challenging to select an appropriate model for characterizing perfusion in specific conditions. Furthermore, an analytical framework that could transform the stretched exponential fit parameters to the intuitive picture of the two-compartment IVIM model could provide utility in comparing and interpreting experimental results.

This work establishes an analytical framework for direct comparisons between the biexponential and stretched exponential models. The explicit representation of the equivalent biexponential model parameters is provided in terms of the stretched exponential model for facilitating comparisons. These analytical expressions and numerical analyses are used to compare models under a range of model parameters that are experimentally encountered in the application of DW-MRI to perfusion measurements of skeletal muscle. The analysis of in vivo DW-MRI data from skeletal muscle is presented as they relate to analytic and numerical results.

2. Definitions and Properties

The most commonly implemented DW-MRI sequence is the pulsed gradient spin echo (PGSE) [13]. This sequence applies diffusion-sensitizing gradients before and after a refocusing radiofrequency pulse separated by a fixed diffusion waiting period, where adjustments to gradient strength are used to set the desired level of diffusion sensitization. Thus, hydrogen nuclei exhibiting motions along the direction of the diffusion sensitizing gradients lose phase coherence proportional to the amount of motion. For the classic Stejskal–Tanner pulse sequence [14], the model used to describe the magnitude of signal attenuation can be expressed in terms of the gradient duration (δ), strength (g), and the time interval between gradient pulses (Δ):

$$y\left(g,\delta ,\mathsf{\Delta }\right)=y_0{e}^{-D{\left(\gamma g\delta \right)}^2\left(\Delta -\frac{\delta }{3}\right)},$$

(1)

where $y_0$ is the equilibrium magnetization in the absence of diffusion sensitizing gradients, D is the apparent diffusion coefficient with units of mm²/s, γ is the gyromagnetic ratio, and Δ$\ge \delta$ [15]. With the majority of clinical MRI scanners using fixed values of δ and Δ, this expression is routinely simplified:

$$y\left(b\right)=y_0{e}^{-Db},$$

(2)

where b is the diffusion-weighting value with units of s/mm². Thus, a common way to adjust diffusion weighting is to vary the gradient amplitude.

The diffusion model has been extended to quantify incoherent motions of blood water protons in the microvascular space. Because the rate of incoherent motions of the blood pool can be orders of magnitude greater than the rate of thermal diffusion of body water, perfusion measurements rely on small values of b, typically sampled in a range between 0 and 200 s/mm². The most well-known perfusion model—IVIM—extends Equation (2) above by taking a biexponential form [6]. This model allows for the fractional pool size and rate of the more rapidly moving blood pool to be quantified separately from that of the tissue water in the extra-vascular space. A more recent emerging perfusion model—FFD--takes the form of a stretched exponential signal model. This model imposes a fractional time derivative on Fick’s law which allows for the modeling of the non-Gaussian jump length statistics expected from blood perfusion within the multiscale vascular bed [16].

2.1. Biexponential Model

As introduced above, the IVIM model uses a biexponential signal model:

$$y_{be}\left(b\right)=y_0\left[f_p{e}^{-D_{p}b}+\left(1-f_p\right)e^{-Db}\right],$$

(3)

where $y_{be}\left(b\right)$ is the MR signal intensity at a particular $b$ value, $b$ is the diffusion-sensitizing factor with units $\mathrm{s}/{\mathrm{mm}}^2$, $y_0$ is the equilibrium magnetization when $b=0$, $f_p$ is the perfusion fraction ($0\le f_p\le 1$), $D_p$ is the pseudodiffusion coefficient (or perfusion coefficient) with units ${\mathrm{mm}}^2/\mathrm{s}$, and $D$ is the apparent diffusion coefficient of water molecules in the extra-vascular tissue with units ${\mathrm{mm}}^2/\mathrm{s}$. The two exponentials in the model correspond to compartments of intra- and extra-vascular tissue that are characterized by the decay rates, $D_p$ and $D$, respectively. Extending the number of compartments with an improved multi-exponential model can account for additional compartments in the vascular tissue, but it is challenged by the increased number of fitting parameters.

2.2. Stretched Exponential Model

The stretched exponential signal model is a parsimonious representation of anomalous superdiffusion, which is a non-Gaussian diffusion that exhibits greater particle spread than would be expected from Gaussian diffusion [17]. This signal model is expressed as a function of the $b$ value as follows:

$$y_{se}\left(b\right)=y_0{e}^{(-{\left(D_{se}·b\right)}^{\alpha /2})}$$

(4)

where $D_{se}$ is the diffusion coefficient in the stretched exponential model with units ${\mathrm{mm}}^2/\mathrm{s}$, and $\alpha$ is the stretching exponent representing the fractional derivative order [17]. Note that, to exhibit sensitivity to perfusion in a complex and heterogeneous system like skeletal muscle, the stretched exponential model requires the sampling of b so that $b·D_{se}\le 1$. For diffusion sequence parameters in the range of $b·D_{se}>1$, this model is interpreted in the context of hindered diffusion which probes a different aspect of the tissue [18]. In the derivation of the FFD model, $0<\alpha \le 2$ with smaller values of $\alpha$ reflecting an increase in perfusion [16]. In particular, the stretched exponential becomes a monoexponential model when $\alpha =2$. In the remainder of the current work, we use a simplifying step using $\widehat{\alpha }=\alpha /2$ where $0<\widehat{\alpha }\le 1$.

Although it is not the physical picture of the FFD model, another interpretation of the stretched exponential model is that it represents a sum of exponentials [19,20]. Thus, numerically the biexponential model can be considered as a truncated representation of the stretched exponential model. Compared with the biexponential model, the parsimonious stretched exponential model has one less fitting parameter, leading to an inherently smaller fitting estimate variance [21].

3. Materials and Methods

3.1. Biexponential Representation of the Stretched Exponential Model Parameters

To obtain the biexponential representation of the stretched exponential model, we consider the three biexponential model parameters $\left\{f_p, D_p,D\right\}$ as unknowns, with $y_0$ being eliminated from Equation (3) by normalizing diffusion weighted images by the b = 0 image. The goal is to solve for these three unknowns given the two stretched exponential parameters $\left\{\widehat{\alpha },D_{se}\right\}$. Three unknowns require three samples from the stretched exponential signal to form a fully determined case. When there are more than three samples, an over-determined case is established, leading to an estimate of biexponential parameter values rather than exact analytical solutions. For both models, $b$ values are sampled either uniformly spaced or non-uniformly spaced.

The stretched exponential model in Equation (4) can be normalized by the equilibrium magnetization at the $b=0$ condition, $y_0$:

$${\tilde{y}}_{se}\left(b\right)=\frac{y\left(b\right)}{y_0}=e^{\left(-{\left(D_{se}·b\right)}^{\widehat{\alpha }}\right)}={(e^{-D_{se}^{\alpha }})}^{b^{\widehat{\alpha }}}$$

(5)

Similarly, the biexponential model in Equation (3) can be normalized to:

$${\tilde{y}}_{be}\left(b\right)=\frac{y\left(b\right)}{y_0}=f_p{e}^{-D_{p}b}+\left(1-f_p\right)e^{-Db}=f_p{e}^{-D_{p}b}+f{e}^{-Db}$$

(6)

with $f=\left(1-f_p\right)$ and $b\in \left\{0,1,2,\dots \right\}$. Given the sampling interval ($\tau$) of $b$, $D_p={\widehat{D}}_p\tau$ and $D=\widehat{D}\tau$ where ${\widehat{D}}_p$ and $\widehat{D}$ are the respective pseudodiffusion coefficient and apparent diffusion coefficient normalized by $\tau$.

Making the substitution $\eta =e^{-D_{se}^{\widehat{\alpha }}}$ and taking the z-transform of ${\tilde{y}}_{se}\left(b\right)$ in Equation (5) and that of ${\tilde{y}}_{be}\left(b\right)$ in Equation (6) leads to:

$$Y_{se}\left(z\right)={\displaystyle \sum }_{b=0}^{\infty }{\eta }^{b^{\widehat{\alpha }}}{z}^{-b}=1+\eta z^{-1}+{\eta }^{2^{\widehat{\alpha }}}{z}^{-2}+{\eta }^{3^{\widehat{\alpha }}}{z}^{-3}+\dots$$

(7)

and

$$Y_{be}\left(z\right)=\frac{f_p}{1-e^{-D_p}{z}^{-1}}+\frac{f}{1-e^{-D}{z}^{-1}}$$

(8)

Expanding and rewriting Equation (8) using ${\rho }_1=e^{-D_p}$ and ${\rho }_2=e^{-D}$ yields:

$$Y_{be}\left(z\right)=\frac{b_0+b_1{z}^{-1}}{1+a_1{z}^{-1}+a_2{z}^{-2}}$$

(9)

where

$$b_0=f_p+f=1;b_1=-\left({\rho }_{1}f+{\rho }_2{f}_p\right);a_1=-\left({\rho }_1+{\rho }_2\right);\mathrm{and}{a}_2={\rho }_1{\rho }_2.$$

(10)

Long division to the right-hand-side (RHS) of Equation (9) results in:

$$\begin{array}{c}{Y}_{be}\left(z\right)=b_0+\left(b_1-b_0{a}_1\right)z^{-1}-\left\{b_0{a}_2+a_1\left(b_1-b_0{a}_1\right)\right\}{z}^{-2}\\ +\left[a_1\left\{b_0{a}_2+a_1\left(b_1-b_0{a}_1\right)\right\}-a_2\left(b_1-b_0{a}_1\right)\right]z^{-3}+\dots \end{array}$$

(11)

Now the RHS of Equations (11) and (7) are both in the form of the $z$-term sums. Comparing coefficients of corresponding z-terms of Equations (11) and (7), we obtain:

$$\begin{array}{c}{b}_0=1,\\ \left(b_1-b_0{a}_1\right)=\eta ,\\ -\left\{b_0{a}_2+a_1\left(b_1-b_0{a}_1\right)\right\}={\eta }^{2^{\widehat{\alpha }}},\\ a_1\left\{b_0{a}_2+a_1\left(b_1-b_0{a}_1\right)\right\}-a_2\left(b_1-b_0{a}_1\right)={\eta }^{3^{\widehat{\alpha }}},\\ ⋮\end{array}$$

(12)

The unknowns $\left\{b_1,a_1,a_2\right\}$ in Equation (12) are associated with the biexponential model parameters $\left\{f_p, D_p,D\right\}$ according to Equation (10), and can be solved with a minimum of three stretched exponential model samples $\left\{\eta , {\eta }^{2\widehat{\alpha }},{\eta }^{3\widehat{\alpha }}\right\}$, as displayed in Equation (12). As implied from Equation (12), these analytical expressions based on the z transform require uniform $b$ value sampling.

In the case where three uniformly spaced samples of the stretched exponential model are given, the fully determined case is established. This yields exact solutions to $b_1, a_1,$ and $a_2$ and, subsequently, exact analytical solutions to $f_p, D_p$ and $D$. By simplifying the Equation (12) by back-substituting terms, we obtain:

$$\{\begin{array}{c}-a_2-a_1\eta ={\eta }^2{}^{\widehat{\alpha }},\\ -a_2\eta -a_1{\eta }^{2^{\alpha }}={\eta }^{3^{\widehat{\alpha }}},\end{array}$$

(13)

which can be written in matrix form as:

$$\left[\begin{array}{cc}-\eta & -1\\ -{\eta }^{2^{\alpha }}& -\eta \end{array}\right]\left[\begin{array}{c}{a}_1\\ a_2\end{array}\right]=\left[\begin{array}{c}{\eta }^{2^{\widehat{\alpha }}}\\ {\eta }^{3^{\widehat{\alpha }}}\end{array}\right]$$

(14)

The fully determined 2 × 2 matrix in Equation (14) allows for analytical solutions to $a_1,a_2$, using $\left(b_1-b_0{a}_1\right)=\eta$ in Equation (12), that of $b_1$:

$$a_1=-\frac{{\eta }^{3^{\widehat{\alpha }}}-\eta {\eta }^{2^{\widehat{\alpha }}}}{{\eta }^{2^{\alpha }}-{\eta }^2};a_2=\frac{\eta {\eta }^{3^{\widehat{\alpha }}}-{\eta }^{2^{\widehat{\alpha }+1}}}{{\eta }^{2^{\alpha }}-{\eta }^2};b_1=\frac{\eta \left(2{\eta }^{2^{\widehat{\alpha }}}-{\eta }^2\right)-{\eta }^{3^{\widehat{\alpha }}}}{{\eta }^{2^{\widehat{\alpha }}}-{\eta }^2}.$$

(15)

Then, $a_1$ and $a_2$ from Equation (15) are substituted in Equation (10) to solve for ${\rho }_1$ and ${\rho }_2$:

$${\rho }_{1,2}=\frac{-a_1\pm \sqrt{a_1^2-4a_2}}{2}$$

(16)

For the IVIM model, ${\rho }_1=e^{-D_p}$, ${\rho }_2=e^{-D}$, and the pseudo-diffusion coefficient $D_p$ is larger than the diffusion coefficient $D$, thus, ${\rho }_1<{\rho }_2$. We obtain:

$$\{\begin{array}{c}{\rho }_1=\frac{-a_1-\sqrt{a_1^2-4a_2}}{2},\\ {\rho }_2=\frac{-a_1+\sqrt{a_1^2-4a_2}}{2},\end{array}$$

(17)

The values $f_p$ and $f$ can be determined based on Equation (10) where:

$$\{\begin{array}{c}{b}_0=1=f_p+f,\\ b_1=-\left({\rho }_{1}f+{\rho }_2{f}_p\right).\end{array}$$

(18)

This appears as follows in matrix form:

$$\left[\begin{array}{cc}1& 1\\ {\rho }_2& {\rho }_1\end{array}\right]\left[\begin{array}{c}{f}_p\\ f\end{array}\right]=\left[\begin{array}{c}1\\ -b_1\end{array}\right]$$

(19)

with the following solutions to $f_p$ and $f$:

$$\{\begin{array}{c}{f}_p=\frac{b_1+{\rho }_1}{{\rho }_1-{\rho }_2},\\ f=1-f_p=-\frac{b_1+{\rho }_2}{{\rho }_1-{\rho }_2}.\end{array}$$

(20)

Our goal is to calculate the biexponential model parameters $\{f_p$, $D_p,D\}$ based on the known stretched exponential parameters $\left\{\widehat{\alpha ,}{D}_{se}\right\}$. Up until this point, the following relationships between stretched and biexponential model parameters have been obtained:

$$\{\begin{array}{c}{f}_p=\frac{b_1+{\rho }_1}{{\rho }_1-{\rho }_2}\\ e^{-D_p}={\rho }_1=\frac{-a_1-\sqrt{a_1^2-4a_2}}{2}\\ e^{-D}={\rho }_2=\frac{-a_1+\sqrt{a_1^2-4a_2}}{2}\end{array}$$

(21)

Combining solutions to $b_1$, ${\rho }_1$, ${\rho }_2$, $a_1$ and $a_2$ into Equation (21), produces analytical expressions for $\{f_p$, $D_p,D\}$ in the form of $\eta$ which is defined above as $\eta =e^{-D_{se}^{\widehat{\alpha }}}$.

$$\begin{array}{c}{D}_p=-\ln \{\frac{{\eta }^{3^{\widehat{\alpha }}}-\eta {\eta }^{2^{\widehat{\alpha }}}-\sqrt{{\eta }^{2*3^{\widehat{\alpha }}}+\eta {\eta }^{3^{\widehat{\alpha }}}\left(4{\eta }^2-6{\eta }^{2^{\widehat{\alpha }}}\right)+4{\eta }^{3*2^{\widehat{\alpha }}}-3{\eta }^2{\eta }^{2^{\widehat{\alpha }+1}}}}{2\left({\eta }^{2^{\widehat{\alpha }}}-{\eta }^2\right)}\},\\ D=-\ln \left\{\frac{{\eta }^{3^{\widehat{\alpha }}}-\eta {\eta }^{2^{\widehat{\alpha }}}+\sqrt{{\eta }^{2*3^{\widehat{\alpha }}}+\eta {\eta }^{3^{\widehat{\alpha }}}\left(4{\eta }^2-6{\eta }^{2^{\widehat{\alpha }}}\right)+4{\eta }^{3*2^{\widehat{\alpha }}}-3{\eta }^2{\eta }^{2^{\widehat{\alpha }+1}}}}{2\left({\eta }^{2^{\widehat{\alpha }}}-{\eta }^2\right)}\right\},\\ f_p=0.5\{1+\frac{{\eta }^{3^{\widehat{\alpha }}}-\eta \left(3{\eta }^{2^{\widehat{\alpha }}}-2{\eta }^2\right)}{\sqrt{{\eta }^{2*3^{\widehat{\alpha }}}+\eta {\eta }^{3^{\widehat{\alpha }}}\left(4{\eta }^2-6{\eta }^{2^{\widehat{\alpha }}}\right)+4{\eta }^{3*2^{\widehat{\alpha }}}-3{\eta }^2{\eta }^{2^{\widehat{\alpha }+1}}}}\},\\ f=1-f_p=0.5\{1-\frac{{\eta }^{3^{\widehat{\alpha }}}-\eta \left(3{\eta }^{2^{\widehat{\alpha }}}-2{\eta }^2\right)}{\sqrt{{\eta }^{2*3^{\widehat{\alpha }}}+\eta {\eta }^{3^{\widehat{\alpha }}}\left(4{\eta }^2-6{\eta }^{2^{\widehat{\alpha }}}\right)+4{\eta }^{3*2^{\widehat{\alpha }}}-3{\eta }^2{\eta }^{2^{\widehat{\alpha }+1}}}}\}.\end{array}$$

(22)

3.2. Biexponential Estimates of the Stretched Exponential Model Parameters

In practice, DW-MRI images are acquired using several non-uniformly spaced $b$ values (i.e., >3) for estimating perfusion and diffusion properties. This sampling scheme results in an over-determined matrix for solving $a_1$, $a_2$, $b_1$, and subsequently $f_p, D_p$ and $D$. Specifically, the matrix in the left-hand side of Equation (14) becomes:

$$\left[\begin{array}{c}\begin{array}{c}\begin{array}{cc}-{\eta }^{b_1^{\widehat{\alpha }}}& -{\eta }^{b_0^{\widehat{\alpha }}}\\ -{\eta }^{b_2^{\widehat{\alpha }}}& -{\eta }^{b_1^{\widehat{\alpha }}}\end{array}\\ ⁝\end{array}\\ \begin{array}{cc}-{\eta }^{b_{k-1}^{\widehat{\alpha }}}& -{\eta }^{b_{k-2}^{\widehat{\alpha }}}\end{array}\end{array}\right]\left[\begin{array}{c}{a}_1\\ a_2\end{array}\right]=\left[\begin{array}{c}\begin{array}{c}{\eta }^{b_2^{\widehat{\alpha }}}\\ {\eta }^{b_3^{\widehat{\alpha }}}\end{array}\\ ⁝\\ {\eta }^{b_k^{\widehat{\alpha }}}\end{array}\right]$$

(23)

where $b_k=k$ for $k\in \left\{0,1,2,\dots \right\}$ in a generalized form. Equation (23) can be written as $\Pi a=p$, where $\Pi$ is a $\left(k-1\right)\times 2$ tall matrix, $a$ is a column vector consists of $a_1$ and $a_2$, and $p$ is a column vector with $\left(k-1\right)$ elements. The method of least squares (LS) [22] to obtain $\widehat{a}={\left({\Pi }^{\mathrm{T}}\Pi \right)}^{-1}{\Pi }^{\mathrm{T}}p$, as estimates of $a$. The analytical expression for the biexponential fit estimates that $\widehat{f_p}$, ${\widehat{D}}_p$ and $\widehat{D}$, can be computed in a manner similar to subsequent steps following Equation (15) in Section 3.1, though the expressions may not be as tractable as Equation (22), obtained for the fully determined case.

3.2.1. Error Analysis

Mean squared error (MSE) in the biexponential modeling of the stretched exponential process can be expressed as:

$$E\left\{e^2\left(b\right)\right\}=E{\left\{{\tilde{y}}_{se}\left(b\right)-{\tilde{y}}_{be}\left(b\right)\right\}}^2$$

(24)

where $e\left(b\right)$ is the model representation error, and $E$ is the expectation operator. Expanding the right-hand-side of (24) leads to:

$$E\left\{e^2\left(b\right)\right\}=E\left\{{\tilde{y}}_{se}^2\left(b\right)+{\tilde{y}}_{be}^2\left(b\right)-2{\tilde{y}}_{se}\left(b\right){\tilde{y}}_{be}\left(b\right)\right\}.$$

(25)

Substituting respective values of ${\tilde{y}}_{se}\left(b\right)$ and ${\tilde{y}}_{be}\left(b\right)$ from Equations (5) and (6) into (25) leads to:

$$E\left\{e^2\left(b\right)\right\}=E\left\{{\eta }^{2b^{\widehat{\alpha } }}\right\}+E\left\{{\left(f_p{\rho }_1^b+f{\rho }_2^b\right)}^2\right\}-2E\left\{{\eta }^{b^{\widehat{\alpha } }}\left(f_p{\rho }_1^b+f{\rho }_2^b\right)\right\}$$

(26)

with $\eta =e^{-D_{se}^{\widehat{\alpha }}}$, and ${\rho }_1=e^{-D_p}$ and ${\rho }_2=e^{-D}$. Equation (26) shows that MSE is a function of the stretching exponent $\widehat{\alpha }$. MSE is zero when $f_p{\rho }_1^b={\eta }^{b^{\widehat{\alpha }}}$ and $f{\rho }_2^b=0$, or vice versa, implying $\widehat{\alpha }=0,$1. The stretched exponential becomes a constant value signal for $\widehat{\alpha }=0$, and a mono-exponential signal for $\widehat{\alpha }=1$.

3.2.2. Error Analysis in the Presence of Noise

In the presence of the independent and identically distributed additive Gaussian noise $v\left(b\right)$ of zero mean and ${\sigma }^2$ variance, Equation (26) becomes:

$$E\left\{e^2\left(b\right)\right\}=E\left\{{\eta }^{2b^{\widehat{\alpha } }}\right\}+E\left\{{\left(f_p{\rho }_1^b+f{\rho }_2^b\right)}^2\right\}-2E\left\{{\eta }^{b^{\widehat{\alpha } }}\left(f_p{\rho }_1^b+f{\rho }_2^b\right)\right\}+{\sigma }^2$$

(27)

with $E\left\{v^2\left(b\right)\right\}={\sigma }^2$. MSE is reduced to ${\sigma }^2$ when $\widehat{\alpha }=0,$1. This will serve as the lower bound for the MSE. For $0<\widehat{\alpha }<1$, MSE will be larger than ${\sigma }^2$.

3.2.3. Model Simulations and Cross-Model Error Analysis

To compare the analytical and numerical solutions derived in Section 3.1 and Section 3.2, forward simulation of stretched exponential signal decays are generated and fitted using the biexponential expressions. Noise is not considered for either the forward or estimated model. All simulations are implemented with custom-written MATLAB (MathWorks 2019b, Natick, MA, USA) scripts.

The z-plane analysis is performed as a way of examining the relationship between biexponential parameters to the stretched exponential signal model. A simulation is implemented for the fully determined case when $b$ = 0, 1, 2, and 3 and $D_{se}=1/b_{max}$. This utilizes the analytical expressions relating the biexponential model parameters to that of the stretched exponential, based on the fully determined case. These $b$ values, although not representative of values used experimentally, in combination with the selection of $D_{se}=1/b_{max}$, reflect the typical scaling used in experimental settings $D_{se}·b_{max}=1$. With this analysis, analytical expressions obtained for biexponential model parameters $\left\{f_p, f, D_p, D\right\}$ in Equation (22) can be jointly visualized as a function of stretched exponential parameters $\left\{\widehat{\alpha }, D_{se}\right\}$, using the so-called three-dimensional z-plane plot. According to Equation (8), the biexponential model corresponds to two poles, ${\rho }_1=e^{-D_p}$ and ${\rho }_2=e^{-D}$ with their respective damping amplitudes $f_p$ and $f$. In the three-dimensional z plane plot, the two in-plane axes refer to the real and imaginary parts of a pole, and the vertical axis refers to the corresponding damping amplitude. When $D_{se}$ is fixed, the trajectory of a pole is directed by the stretched exponential model parameter $\widehat{\alpha }$.

Simulations for the over-determined case compare estimated biexponential parameters with their corresponding analytical solutions for $0<\widehat{\alpha }\le 1$, at a fixed $b_{max}$ and $D_{se}=\frac{1}{b_{max}}$. Sampling schemes of $b$ values at various sample spacings were examined and compared with the analytical solutions to biexponential model parameters. Analytical expressions are computed and shown for two limiting cases: (1) $D_{se}=1/3$ and (2) $D_{se}=1/400$.

The divergence between the simulated signal $y_{se}\left(b\right)$ and estimated signal ${\widehat{y}}_{se}\left(b\right)$ of the stretched exponential model is quantified using mean-squared-error (MSE) analysis. Given uniformly spaced b values and $D_{se}=\frac{1}{b_{max}}$, MSE is computed as a function of $\widehat{\alpha }$, defined as $\epsilon \left(\widehat{\alpha }\right)=\frac{1}{N}{\displaystyle \sum }_{n=1}^N{\left(y_{se}\left(b_n\right)-{\widehat{y}}_{se}\left(b_n\right)\right)}^2$, where $N$ is the number of b values, and signal values $y_{se}\left(b_n\right)$ and ${\widehat{y}}_{se}\left(b_n\right)$ are modulated by the argument $\widehat{\alpha }$ (refer to equation (5)). Additionally, the autocorrelation of the error sequence $\epsilon \left(\widehat{\alpha },m\right)=\frac{1}{N}{r}_{ee}\left(m\right)$, is computed as is a function of $\widehat{\alpha }$ and delay $m$, where $m$ is the b-value lag, and $r_{ee}\left(m\right)$ represents the autocorrelation function of the residual $e\left(b_n\right)=y_{se}\left(b_n\right)-{\widehat{y}}_{se}\left(b_n\right)$.

3.3. In Vivo MRI Acquisition

To illustrate the analytical and numerical analysis of the IVIM and FFD models, we compared model fit parameters from DW-MRI data from a single healthy volunteer. Data were obtained as part of the Genetic and Epigenetic Signatures of Translational Aging Laboratory Testing (GESTALT) study under approval by IRB through the National Institute of Environmental Health Sciences.

The exercise protocol was described in existing publications [16,23]. Briefly, DW-MRI was obtained before and immediately after performing an in-magnet knee extension exercise with an MR-compatible ergometer (Quadspect; Ergospect, Innsbruck, Austria) at a frequency of 45 repetitions per minute for 2.5 min using the left leg. The pedal resistance was fixed at 0.4 bar.

The participant was imaged using a 3T Achieva MRI scanner (Phillips Healthcare, Best, Netherlands) using a two-channel SENSE flexible coil positioned at the mid-thigh level. The DWI measurements were acquired using a single transverse slice of thickness 22 mm at mid-thigh, using a spin-echo single-shot EPI readout with diffusion sensitization in the slice direction, with triple fat suppression to reduce the effect of fat signal on the measurements. Sequence parameters included TR = 3000 ms, TE = 46 ms, SENSE factor = 2, FOV = 256 × 225 mm, matrix size = 224 × 224, partial Fourier factor = 0.5, diffusion gradient duration δ = 27 ms, diffusion gradient interval Δ = 35 ms, and 14 non-uniformly spaced $b$ values: 0, 3, 7, 10, 15, 20, 25, 30, 40, 50, 70, 100, 200, and 400 s/mm². DWI acquisition at rest performed 20 signal averages (acquisition time = 16.05 min) in order to increase SNR under low perfusion conditions, while post-exercise acquisition performed 4 signal averages (acquisition time = 3.25 min). All images exhibited SNR > 200, deemed sufficient for robust estimates of parameters [23,24].

Regions of interest (ROI) were selected within the rectus femoris and vastus lateralis, which were part of the active muscle groups during knee extension exercise. ROIs were drawn to exclude blood vessels, fat and fasciae. Pixel-wise fits were performed within each ROI using the biexponential (Equation (6)) and stretched exponential (Equation (5)) models, respectively. The connection between parameters from the two model fits was compared with the simulation results obtained from the overdetermined case described in Section 3.2. Parameter comparisons were made from resting muscle and muscle after exercise, representing a large dynamic range of model parameters observed in vivo.

4. Results

4.1. Fully-Determined Case

Analytical expressions relating the biexponential and stretched exponential model parameters were used to generate Figure 1 and Figure 2, using uniformly spaced $b$ values ($b=0$, 1, 2, 3). These analytical relationships for the parameters $f_p, D_p$ and $D$ showed excellent correspondence with the numerical solutions derived in Section 3.2 over the range of simulated $\widehat{\alpha }$ and $D_{se}$ values, respectively.

Figure 3 shows each biexponential model parameter with respect to the bivariant combinations of the stretched exponential model parameters $\widehat{\alpha }$ and $D_{se}$. These plots illustrate how biexponential parameters are influenced by the relationship between $D_{se}$ and $b_{max}$. Specifically, these plots illustrate the relatively more pronounced nonlinear influence of $\widehat{\alpha }$ on biexponential parameters, as $D_{se}·b_{max}$ decreases from 1.

Figure 4 shows the z-plane analysis, which illustrates the pole trajectories and their corresponding damping amplitudes. Trajectories of both poles lie on the real-axis of the z-plane, indicating that the biexponential estimates based on the stretched exponential data reconstruct real exponentials only. Both poles lie within the unit circle, indicating that the biexponential estimates of the stretched exponential model are stable. The trajectories of the two poles do not coincide, implying that the solution space for the biexponential model is localized for estimating the forward model of the stretched exponential model.

In the projection on the z-plane, a pole closer to the unit circle boundary exhibits a slower decay rate that corresponds to a narrower diffusion jump length distribution. The pole closer to the center of the unit circle exhibits a faster decay rate and contributes to the relatively heavy tail in the jump length distribution. As $\widehat{\alpha }$ changes progressively from 0⁺ to 1, the pole associated with pseudo-diffusion travels from the center of the unit circle toward the boundary of the unit circle, while the pole associated with diffusion travels from the boundary toward the center of the unit circle. The amplitude of the pole associated with pseudo-diffusion vanishes, whereas the amplitude of the pole associated with diffusion becomes progressively larger. The pole with initially higher amplitude ($f_p$) accounts for the starting values of the stretched exponential, and the pole with a relatively lower initial amplitude ($f$) accounts for the final values of the stretched exponential.

When $\widehat{\alpha }=1$, the stretched exponential converges to a monoexponential. In this case, the biexponential model becomes redundant, which is reflected on the pole trajectory plot (Figure 4) in the vanishing of the pole that is driven by the pseudo-diffusion process ($f_p=0$). The diffusion signified by the remaining pole is sufficient to characterize the diffusive motion of water molecules.

4.2. Over-Determined Case

Numerical simulations provide a means of comparing the two models under experimentally relevant b-value sampling schemes. Figure 5 shows the biexponential parameters estimated from simulated stretched exponential decays using Prony’s method (described in Section 3.2). These results are obtained using more than four samples from the uniformly spaced DW-MRI signal (including $y_0$ acquired at $b=0$). In Figure 5a,b, numerical solutions progressively deviate from the analytical solutions derived for the fully determined case (Equation (22)) as the spacing of uniformly distributed samples of $b$ values, $\tau$, increases. Interestingly, with the decimation of the b-value sampling, all parameters $\left\{f_p, f,D,D_p\right\}$ assume a near-linear relationship with respect to $\widehat{\alpha }$, particularly at larger values of $\widehat{\alpha }$. The influence of sampling on $D_p$ and $D$ exhibits more complicated relationships with b value decimation. For example, numerical estimates of $D_p$ deviate more with the increased sampling of the b value and approach the analytical result with decimation. Numerical estimates of $D$ show the largest deviations from the analytical expressions taking a different nonlinear shape with large overestimations at larger values of $\widehat{\alpha }$.

4.3. Analysis of Fit Error

For the fully determined case that uses four uniformly spaced $b$ values (including $b=0$) to obtain analytical solutions, zero residual is observed between decay signals from these two models. This result does not indicate that the biexponential model can perfectly represent the stretched exponential, but rather that there is an insufficient number of samples to observe deviations between the stretched exponential decay and the biexponential model.

In the overdetermined case, when more than four uniformly spaced b values (including $b=0$) are used for the simulation, the MSE becomes nonzero, and its value varies with $\widehat{\alpha }$. Figure 6 shows an example of the MSE (panel (a) and the autocorrelation of fit error (panel (b)) over the range of $\widehat{\alpha }$ using the noiseless and overdetermined b-value sampling that satisfies $D_{se}·b_{max}=1$. The MSE (Figure 6a) remains relatively low for larger values of $\widehat{\alpha }$, suggesting difficulty in discriminating between biexponential and stretched exponential models based on goodness of fit. As values of $\widehat{\alpha }$ decrease, the autocorrelation function of error (Figure 6b) shows oscillations with increasing lags, reflecting an increase in the fit error of the biexponential approximation of the stretched exponential function. The cross-section of Figure 6b at lags = 0 corresponds to the all-positive profile of MSE in Figure 6a.

The stretched exponential model, as shown by others, can be represented as a sum of multiple exponentials. Therefore, the biexponential model contains fewer exponentials than the stretched exponential model and for $\widehat{\alpha }<1$ tends to underfit. This can be observed in Figure 6 based on the shape of the MSE as a function of $\widehat{\alpha }$.

4.4. In Vivo Results

Comparisons of model fit parameters derived from experimental data from skeletal muscle demonstrated a strong resemblance to theoretical curves. Specifically, comparisons of biexponential parameters with respect to the stretched exponential parameter $\widehat{\alpha }$ are shown from selected ROI’s in the rectus femoris (Figure 7) and vastus lateralis (Figure 8). The numerical estimates were obtained using the same b values for collecting experiment data and using the mean of $D_{se}$ estimates derived from the stretched exponential fit for the experimental data.

For the perfusion fraction ($f_p$), the experimental data shows a strong correspondence between experimental and theoretical predictions over all $\widehat{\alpha }$ values, which range from 0.2 to 1. Similarly, 1 − f_p shows a strong correspondence between experimental and theoretical values. The pseudo-diffusion coefficient ($D_p$) and diffusion coefficient ($D$) show more scatter around the theoretical curve but largely follow the theoretical prediction. In both muscle groups displayed, these relationships were consistent under physiological conditions at rest and post-exercise hyperemia. Larger deviations of $D_p$ and $D$ in experimental data appear to occur at lower values of $\widehat{\alpha }$. The greater variance in of $D_p$ and $D$ can also be explained by Crammer Rao bound analysis of the biexponential model, which shows that the variance of the diffusion rate constant estimates is directly proportional to their cubed value and the quartic of the sum of $D$ and $D_p$, and inversely related to the quartic of their difference. Thus, the variance in $D_p$ estimates is expected to be large for large rates ($D_p$) and for small separations between $D$ and $D_p$ [25].

5. Discussion and Conclusions

The exact solutions of biexponential parameters in Equation (22) exhibit complicated nonlinear forms in terms of stretched exponential model parameters that are not easily interpreted. However, graphical comparisons provide some key insights that are useful in examining model comparisons, particularly when considering the multi-exponential nature of the stretched exponential function [19]. For example, when assessing only 4 $b$ values spanning the full range of $1/D_{se}$ (Figure 1), the bicomponent estimates appear almost linear with $\widehat{\alpha }$. Under these conditions, both models are indistinguishable based on the assessment of MSE (data not shown), which is not surprising since four data points are insufficient to represent more than two exponentials. Numerical simulations using the increased sampling of $b$ values over the range of $1/D_{se}$ (Figure 5) show how these curves approach more nonlinear relationships between biexponential parameters and $\widehat{\alpha }$. With decreased sample spacing, $f_p$, $1-f_p$ and $D_p$ approach curves are generated by analytical expressions using only the first four b values. In contrast, $D$ still shows a deviation from the analytical expression, which results from the lack of sampling of larger $b$ values in analytical expressions, which limits the amount of information used to approximate slow diffusion rate information. Thus, the analytical expression for the relationship between $D$ and $\widehat{\alpha }$ overfits rapid decay rates (e.g., perfusion information) and poorly represents slow diffusion (e.g., tissue diffusion) information.

Experimental data closely resembles the theoretical curves connecting the biexponential parameters to the stretched exponent. Perfusion ($f_p$) and diffusion ($1-f_p$) fractions from the rectus femoris muscle (Figure 7) fell almost identically on the numerically generated curve. The values $D_p$ and $D$ in the rectus femoris exhibited more noticeable deviations from the numerically generated curve, particularly during hyperemia, where the shape of the curves shows a slightly different slope. This deviation is not as apparent in the vastus lateralis muscle (Figure 8). One possible explanation for this deviation in rectus femoris results is the influence of unsuppressed and chemically shifted fat signal [23]. The single-shot echo-planar MRI readout from this sequence shifts a small amount of unsuppressed fat signal into the rectus femoris muscle creating a bias effect on biexponential fitting [24]. The vastus lateralis muscle, given the orientation of the image field of view and the geometry of the thigh, is not susceptible to this effect. Both $D_p$ and $D$ in the vastus lateralis muscle show more dispersion around the numerical curve over the range of $\widehat{\alpha }$, which is consistent with Crammer Rao lower bound predictions of decay rate estimates in the biexponential model [25,26]. Differences in perfusion profiles between muscle groups both before and after exercise (Figure 7 and Figure 8) are attributed to differences in the demand for oxygenated blood and muscle involvement in knee extension exercise.

Although the focus of this work compares the FFD (i.e., stretched exponential) model with the ubiquitous biexponential form of the IVIM model, it should be noted that an alternative IVIM model was proposed in Le Bihan and colleagues’ original work [6], in which the perfusion term is modeled as an isotopically oriented plug flow. This condition would apply to capillary segments that are longer than the fluid particle step length within the diffusion period. Although outside of the scope of the current work, this assumption could lead to an alternative set of expressions relating the FFD model to the IVIM model. We note that in skeletal muscle, based on the average reported capillary segment length and flow velocity, the estimated time to traverse a capillary segment is ~180–525 ms, suggesting that plasma flow through the capillaries at rest might be better represented using this plug flow assumption [27]. With arteriole flow velocity being greater and segments being longer, it is likely that this kind of model is applicable to microvessels in skeletal muscle and perhaps should be considered in future work.

The problem of model selection is challenging in the space of human DW MRI-derived perfusion and diffusion applications. Models that are physically motivated by underlying physiologic properties (e.g., microvascular blood volume, microvascular blood flow rate, tissue microarchitecture, etc.) are important for leveraging these imaging tools and applying them to important research questions. However, there is a tension between the use of models that precisely represent the physical picture and those that are practical to apply to in vivo data which is inherently limited in, e.g., the sampling of $b$ values and the signal-to-noise ratio. The stretched exponential model has demonstrated much better stability in parameter estimates than the biexponential model [28]. However, its interpretation in the context of physical properties, such as blood pool size, is much less concrete [16]. Thus, quantitative insights about the relationships between parameters from these two models and how these relationships are influenced by the sampling of the diffusion-weighted signal are useful for interpreting these models as applied to experimental data.