Purpose

The Diffusion-Weighted (DW) signal is tipically considered to be a monoexponential decay at low b-values where diffusion follows the statistical mechanics of Brownian motion. However, studies^1,2,3 have shown that at high b-values the DW signal deviates from the monoexponential decay. Consequently, there is significant interest in developing models able to phenomenologically explain this deviation. Multiple models have been proposed to understand the actual in-vivo diffusion process, i.e., the empirical stretched exponential⁴ and the kurtosis⁵ models try to explain these deviations through a heterogeneity and non-gaussianity index. More complex models generalize the Brownian motion, where the Mean Square Displacement (MSD) is proportional to the diffusion time, to the Continuous-Time Random Walk (CTRW) where the MSD is represented as the power law $$$\langle{x}^2(t)\rangle~\approx~t^{\alpha/\beta}$$$ where 𝛼 and 𝛽 are the time and space fractional dimensions of the lattice. This generalization converts Gaussian diffusion (𝛼/𝛽=1) into anomalous diffusion (i.e., super-diffusion when 𝛼/𝛽>1 and subdiffusion when 𝛼/𝛽<1)⁶ depending on the complexity, heterogeneity and viscosity of the tissue⁷. Accordingly, the CTRW model⁷ solves the space-time fractional order diffusion equation, and Magin et. al.⁸ introduced the fractional-time/space order Bloch-Torrey equation for the transverse magnetization ($$$M_{xy}({r},t)$$$):

Eq. 1:$$^C_0D^{\alpha}_t\{M_{xy}({r},t)\}~=~-i\gamma(\vec{{r}}\cdot\vec{{G}}(t))M_{xy}({r},t)~+~D\nabla^{2\beta}M_{xy}({r},t),$$

where $$$^C_0D^{\alpha}_t\{\cdot\}$$$ is the Caputo 𝛼th-fractional order time derivative (0<𝛼<=1), $$$\nabla^{2\beta}$$$ the 2𝛽th-fractional order Laplacian operator in space (1/2<𝛽<=1), $$$\vec{{G}}(t)$$$ the time-varying diffusion-weighting gradients, and D the anomalous diffusion coefficient. However, they only solved the Fractional-Space order Bloch-Torrey (FSBT) equation (𝛼=1, 1/2<𝛽<=1) for Stejskal-Tanner gradients, and the Fractional-Time order Bloch-Torrey (FTBT) equation (0<𝛼<=1, 𝛽=1) for the trivial case of a constant gradient. Therefore, no DW-MRI implementable model has been developed for the FTBT which might better explain the nature of the deviations occurring at high b-values. Thus, in this work we solve the Fractional-Time order Bloch-Torrey equation for generalized smooth diffusion-weighting gradients, and validate the proposed model on synthetic anomalous phantom experiments.

Theory

Given the continuous m^th-order smooth polynomial time window functions (w_m) of Singla et. al.⁹, we constructed diffusion-weighting gradient waveforms of spin-echo DWI as $$$\vec{{G}}(\text{t})~=~g\cdot\text{w}_m(\frac{2t-\delta}{\delta})\vec{{r}}~+~g\cdot\text{w}_m(\frac{2(t-\Delta)-\delta}{\delta})\vec{{r}}$$$, where g, 𝛿, and 𝛥 are the gradient strength, duration, and separation, respectively. Then, we solved the FTBT (Eq. 1 for 0<𝛼<=1, and 𝛽=1) for the case of a single gradient in the z-direction ($$$\vec{{r}}~=~\hat{{z}}$$$ and $$$G_z~=~\vec{{G}}(\text{t})\cdot\hat{{z}}$$$). Firstly, the homogeneous solution is $$$M_{xy}(z,t)~=~M_0\mathrm{E}_\alpha[-i\gamma{z}\cdot{I}^{\alpha}_t\{G_z(t)\}]$$$, where $$$M_{xy}(r,0)~=~M_0$$$, $$$I^{\alpha}_t\{\cdot\}$$$ is the fractional integral of order 𝛼, and $$$\mathrm{E}_\alpha[\cdot]$$$ is the single parameter Mittag-Leffler function¹⁰. Then, by applying truncated versions of the Leibniz rule and the Faá di Bruno’s formula¹¹ for the Caputo operator (note that the m^th-order of the time window functions has to be equal or higher to the index of previous truncations) we obtained the following solution of the Fractional-Time order Bloch-Torrey (FTBT) equation for the transverse magnetization:

Eq. 2:$$M_{xy}(z,t)~=~M_0\mathrm{E}_\alpha[-i\gamma{z}\cdot{I}^{\alpha}_t\{G_z(t)\}]e^{\int_0^tH(\tau)d\tau},$$

where $$H(t)~=~\frac{-^C_0D^{\alpha}_t\{\mathrm{E}_\alpha[-i\gamma{z}\cdot{I}^{\alpha}_t\{G_z(t)\}]\}~-~i\gamma{z}G(t)\mathrm{E}_\alpha[-i\gamma{z}\cdot{I}^{\alpha}_t\{G_z(t)\}]~+~D_\alpha\frac{\partial^2\mathrm{E}_\alpha[-i\gamma{z}\cdot{I}^{\alpha}_t\{G_z(t)\}]}{\partial{z}^2}}{\alpha^C_0D^{\alpha}_t\{\mathrm{E}_\alpha[-i\gamma{z}\cdot{I}^{\alpha}_t\{G_z(t)\}]\}}.$$

Hereinafter, we numerically evaluate this equation for time window functions of m^th-order m=2.

Methods

Differently to Brownian motion where the step size of the random walk varies according to a Gaussian distribution, in continuous-time random walks both the waiting time and the jump length can vary according to different fractional distributions. Therefore, we simulated an anomalous diffusion medium through a random walk where the probability density function (PDF) of the jumps depended on the position (x) of each walker, on an asymmetry unitless parameter (𝜆), and on the diffusion parameter D. We used PERT distributions with PDF ($$$f(x)~=~\frac{(x-A)^{a-1}(C-x)^{b-1}}{\text{B}(a,b)(C-A)^{a+b-1}}$$$) (see Figure 1 for details), where $$$a~=~\frac{5B+C-5A}{C-A},~b~=~\frac{5C-A-4B}{C-A}$$$ and B(a,b) is the beta function. To control the time dependence of the MSD of the simulations, the position of each walker was assigned following a stochastic process with a PERT distribution with parameters A, B, and C dependent on x and the diffusion time as shown in Figure 1. Also, note that at x=0 the distribution was always Gaussian.

We simulated various spin-echo DWI acquisitions with 1,000,000 walkers initially at the center of coordinates (x=0) with D=2x10-3mm2/s, 𝜆 in the range of 0-0.5, and diffusion-weighting gradients with 𝛿=5ms, 𝛥=[11,16,21,26]ms, and g=[0,100,200,300,400,500]mT/m, TE = 2𝛥. b-values of each [𝛥,g] pair are shown in Table 1. Estimation of the anomalous diffusion, 𝛼, and/or 𝛽 coefficient was performed with the proposed model (Eq. 2) and models of Table 2 through a non-linear least-squares estimator. Note that some models have the standard values of 𝛼=1 and/or 𝛽=1.

Results

Figure 2 shows the percentage normalized squared error (NSE) of the different models’ estimations, showing a very good fit for the proposed model (FTBT).

Discussion

This preliminary results illustrate the potential of the FTBT model to estimate anomalous diffusion parameters. It achieves lower NSE than any other model. Interestingly, compared to the other models the FTBT takes advantage of DW data from multiple gradient strengths (g) and separations (𝛥), probably due to the non-linear time dependence of the MSD. Nevertheless, detailed in-vivo DWI validation is required to further understand the anomalous diffusion behavior in tissues.

Conclusion

We have developed a new diffusion-weighting signal model solving the Fractional-Time order Block-Torrey (FTBT) equation under the assumption of anomalous diffusion. The proposed model minimizes the NSE compared to state-of-the-art models on synthetic anomalous phantom experiments. This model may have important applications for in-vivo DW-MRI to characterize complex and heterogeneous tissues.

Acknowledgements

The authors acknowledge grant RTI2018-094569-B-I00 from the Ministerio de Ciencia e Innovación of Spain. The main author would also like to acknowledge the Consejería de Educación of Junta de Castilla y León, and the Fondo Social Europeo for his predoctoral grant.