Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC5365033

Formats

Article sections

Authors

Related links

Front Phys. Author manuscript; available in PMC 2017 March 24.

Published in final edited form as:

PMCID: PMC5365033

NIHMSID: NIHMS851561

Carson Ingo,^{1} Yi Sui,^{2} Yufen Chen,^{3} Todd B. Parrish,^{3} Andrew G. Webb,^{1} and Itamar Ronen^{1,}^{*}

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

The publisher's final edited version of this article is available at Front Phys

In this paper, we provide a context for the modeling approaches that have been developed to describe non-Gaussian diffusion behavior, which is ubiquitous in diffusion weighted magnetic resonance imaging of water in biological tissue. Subsequently, we focus on the formalism of the continuous time random walk theory to extract properties of subdiffusion and superdiffusion through novel simplifications of the Mittag-Leffler function. For the case of time-fractional subdiffusion, we compute the kurtosis for the Mittag-Leffler function, which provides both a connection and physical context to the much-used approach of diffusional kurtosis imaging. We provide Monte Carlo simulations to illustrate the concepts of anomalous diffusion as stochastic processes of the random walk. Finally, we demonstrate the clinical utility of the Mittag-Leffler function as a model to describe tissue microstructure through estimations of subdiffusion and kurtosis with diffusion MRI measurements in the brain of a chronic ischemic stroke patient.

In the first measurements of water diffusion in biological tissue using magnetic resonance imaging (MRI) systems, the term “apparent diffusion coefficient” (ADC) was chosen to highlight the fact that, although the free diffusion coefficient of water at body temperature is ~ 3 × 10^{−3} mm^{2}/s, typical values in white matter (WM) and gray matter (GM) regions of interest in the human brain were found to be an order of magnitude smaller ~ 0.6 – 1.0 × 10^{−3} mm^{2}/s due to hindrances imposed on water self-diffusion by the tissue microstructure [1–3]. Furthermore, it was recognized that the diffusion decay signal does not always follow a monoexponential decay as predicted by the Bloch-Torrey equation when diffusion weighted measurements are sampled at different time and length scales [4]. The non-monoexponential behavior suggested a superposition of two Gaussian diffusion populations with slow and fast ADC values which related to the volume fractions of intracellular and extracellular space within the imaging voxel, known as the biexponential model [5, 6]. Although the biexponential model provides a means to accurately fit the diffusion signal, the estimation of intracellular and extracellular volume fractions does not accurately correspond to the actual physical makeup of the underlying microstructure, either *in vivo* or even in a simple ghost erythrocytes model [7, 8].

As a result, more sophisticated geometrical models have been developed to more subtly identify diffusion properties that provide insight about tissue microstructure [9–13]. Some approaches describe the non-monoexponential signal as a shift from Gaussian diffusion to non-Gaussian diffusion. There have been efforts to utilize the exponential function and stretch the argument by a parameter in the exponent to provide an index that classifies the deviation from Gaussianity [14–17]. Additionally, there has been an attempt to characterize fixed biological tissue with higher moment analysis of the diffusion propagator in order to extract fractal-dimension measures [18, 19]. Another way to identify non-Gaussian diffusion, known as diffusional kurtosis imaging (DKI), uses a Taylor-series expansion of the argument in the exponential function in order to estimate excess kurtosis of the measured signal vs. the Gaussian case of a monoexponential decay [20]. However, due to the parabolic form of the fitting function in DKI, a limit must be placed on the maximum diffusion gradient strength for the fitting function to monotonically decrease with increased diffusion weighting [21]. It is possible to extend the gradient limit by including higher order terms in the expansion, however this requires the inclusion of terms greater than the second and fourth order cumulants that correspond to the second and fourth moments [22]. As very high strength gradients have recently become available, the mathematical form in DKI limits the ability to interrogate tissue microstructure [23]. Recently, the continuous time random walk (CTRW) has been applied to neural tissue to connect a mathematical model with a physical interpretation of diffusion [24–28]. As pointed out in Jensen and Helpern [21], the stretched exponential model in Bennett et al. [14] is not compatible with DKI as it is not an analytic function. However, by considering a model for subdiffusive processes, it is possible to compute the excess kurtosis using moment analysis of the Mittag-Leffler function (MLF). In CTRW theory, the MLF is a generalization of the exponential function, which is an analytic and monotonically decreasing function for all arguments. The MLF, while rigorously defined as a convergent power series, is, nevertheless, challenging to compute and fit for clinical MRI data, which is constrained in sampling due to practical limits on scan time. Here, we formulate simplified fitting forms of the MLF, connect subdiffusion to kurtosis, provide simulations of random walk conditions to illustrate the diffusion physics, and demonstrate measurements of non-Gaussian diffusion in the brain of a chronic stroke patient.

The fundamental concept in CTRW theory is to extend the diffusion equation such that the fractional-order partial derivatives can be utilized as mathematical operators to describe the diffusion propagator, *P*(*x*, *t*):

$$\frac{{\partial}^{\alpha}P(x,t)}{\partial {t}^{\alpha}}={D}_{\alpha ,\beta}\frac{{\partial}^{\beta}P(x,t)}{\partial {\mid x\mid}^{\beta}},$$

(1)

where * ^{α}*/

The justification for making use of the fractional derivative operators is to provide a mathematical means to interpolate from homogeneous and relatively simple systems that exhibit local, Gaussian behavior to heterogeneous and relatively complex systems that exhibit non-local, power-law behavior [25, 30–34]. In the CTRW context, the fractional order operators, *α* and *β*, provide a description of a random walk’s likelihood to have broader distributions of waiting times and jump lengths, respectively, in comparison to classical Brownian motion. When *α* = 1 and *β* = 2, Equation (1) simplifies to the integer-order partial differential equation to describe Gaussian diffusion. Through the mean-squared displacement (MSD), when the ratio 2*α*/*β* < 1 the dynamics are sub diffusive and when 2*α*/*β* > 1 the dynamics are super diffusive [24, 25]. In the most general case in which *α* and *β* are of arbitrary orders, Equation (1) can be readily transformed to Fourier-Laplace space and represented in closed form as the MLF [32, 35]. The Fourier-Laplace transform, *P*(*x*, *t*)→ *p*(*q*, *s*), in Equation (1) is,

$$p(q,s)=\frac{1}{s+{D}_{\alpha ,\beta}{s}^{1-\alpha}{\mid q\mid}^{\beta}}.$$

(2)

Applying the inverse Laplace transformation to Equation (2), we obtain a simple expression for the characteristic function (CF) of the diffusion propagator as,

$$p(q,t)={E}_{\alpha}[-{D}_{\alpha ,\beta}{\mid q\mid}^{\beta}{t}^{\alpha}],$$

(3)

where *E _{α}* is the single-parameter MLF [36, 37]. When

$$p(q,t)=exp(-{D}_{1,2}{\mid q\mid}^{2}t),$$

(4)

The MLF is an entire function defined as a power series expansion,

$$f(z)={E}_{\alpha}(z)=\sum _{k=0}^{\infty}\frac{{z}^{k}}{\mathrm{\Gamma}(\alpha k+1)},$$

(5)

where the Γ function is the generalized form of the factorial function, defined for real numbers [38]. When *α* = 1, through the identity Γ (*k* + 1) = *k* !, Equation (5) becomes,

$$f(z)={E}_{1}(z)=\sum _{k=0}^{\infty}\frac{{z}^{k}}{k!},$$

(6)

which is the Taylor series definition of the exponential function.

For the case of time-fractional subdiffusion (i.e., 0 < *α* ≤ 1 and *β* = 2), in which the distribution of waiting times for the random walk follows power-law behavior (that is, more likely to wait longer times prior to each step), but the distribution of step lengths is Gaussian, Equation (2) becomes,

$$p(q,s)=\frac{1}{s+{D}_{\alpha ,2}{s}^{1-\alpha}{\mid q\mid}^{2}},$$

(7)

and applying the inverse Laplace transform yields,

$$p(q,t)={E}_{\alpha}[-{D}_{\alpha ,2}{\mid q\mid}^{2}{t}^{\alpha}].$$

(8)

For the case of space-fractional superdiffusion (i.e., *α* = 1 and 1 < *β* ≤ 2), in which the distribution of waiting times is Gaussian, but the distribution for the random walker’s step lengths follows power-law behavior (that is, more likely to jump further), Equation (2) becomes,

$$p(q,s)=\frac{1}{s+{D}_{1,\beta}{\mid q\mid}^{\beta}},$$

(9)

and applying the inverse Laplace transform yields,

$$p(q,t)=exp(-{D}_{1,\beta}{\mid q\mid}^{\beta}t).$$

(10)

In order to write *D _{α}*

$${D}_{\alpha ,\beta}\equiv {D}_{1,2}\frac{{\tau}^{1-\alpha}}{{\mu}^{2-\beta}},$$

(11)

where μ (e.g., μm) and *τ* (e.g., ms) are heuristic parameters that can be estimated to preserve the units for the diffusion coefficient, however, similar parameters have been for conservation of mass and heavy tailed limit convergence in fractal and fractional dynamics [33, 39–41]. When *α* = 1 and *β* = 2, Equation (11) shows that *D _{α}*

$$b\equiv {(\gamma G\delta )}^{2}t,$$

(12)

where *γ* is the gyromagnetic ratio, *G* is the amplitude of the diffusion gradient, *δ* is the duration of the diffusion gradient pulse, and *t* is the effective diffusion time. The parameter *q* is defined as,

$$q\equiv \gamma G\delta .$$

(13)

For the case of the spin-echo variant of the Stejskal-Tanner diffusion-weighted pulse sequence,

$$t\equiv \mathrm{\Delta}-\frac{\delta}{3},$$

(14)

where Δ is the time between diffusion gradient pulses, such that,

$$b\equiv {(\gamma G\delta )}^{2}(\mathrm{\Delta}-\frac{\delta}{3})={q}^{2}t.$$

(15)

Combining Equations (3) and (11),

$$p(q,t)={E}_{\alpha}[-{D}_{1,2}\frac{{\tau}^{1-\alpha}}{{\mu}^{2-\beta}}{q}^{\beta}{t}^{\alpha}],$$

(16)

which can also be written as,

$$p(q,t)={E}_{\alpha}[-{D}_{1,2}\frac{{\tau}^{1-\alpha}}{{\mu}^{2-\beta}}{({q}^{2})}^{\beta /2}{t}^{\alpha}].$$

(17)

By substituting *q*^{2} with *b*/*t*, Equation (17) becomes,

$$p(\frac{b}{t},t)={E}_{\alpha}[-{D}_{1,2}\frac{{\tau}^{1-\alpha}}{{\mu}^{2(1-\beta /2)}}{(\frac{b}{t})}^{\beta /2}{t}^{\alpha}],$$

(18)

which can be rearranged as,

$$p(\frac{b}{t},t)={E}_{\alpha}[-{D}_{1,2}\frac{{\tau}^{1-\alpha}}{{\mu}^{2(1-\beta /2)}}{(b)}^{\beta /2}{t}^{(\alpha -\beta /2)}],$$

(19)

and the substitution can be made to define the classical diffusion coefficient, *D*, raised to the order *β*/2 as,

$${D}^{\beta /2}\equiv {D}_{1,2}\frac{{\tau}^{1-\alpha}}{{\mu}^{2(1-\beta /2)}}{t}^{(\alpha -\beta /2)},$$

(20)

in which *D ^{β}*

$$p(b)={E}_{\alpha}[-{(bD)}^{\beta /2}],$$

(21)

such that *D*, *α*, and *β* can be estimated from a minimum of three data samples. Equation (21) is a compact equation with which to interrogate diffusion and perfusion properties in biological tissues in the realms of sub-, super-, and Gaussian diffusion to capture, for example, intravascular incoherent motion (IVIM) at small *b*-values [42]. However, if one is interested in capturing only subdiffusion and Gaussian diffusion properties of WM and GM of the brain, which was demonstrated to be the general behavior in Ingo et al. [24] as *β* → 2 and 2*α*/*β* ≤ 1 for most voxels, Equation (21) can be further condensed by setting *β* = 2,

$$p(b)={E}_{\alpha}(-bD).$$

(22)

As an alternative approach to develop a model for subdiffusion, Equations (8) and (11) can be combined by setting *β* = 2 to give,

$$p(q,t)={E}_{\alpha}[-{D}_{1,2}{\tau}^{1-\alpha}{\mid q\mid}^{2}{t}^{\alpha}],$$

(23)

which for estimated values of *τ* *t*, as demonstrated in Ingo et al. [24], produces an effective mathematical form given by the two parameter model in Equation (22) such that the equations are interchangeable under the condition *τ* *t*.

Qualitatively, kurtosis is a measure to non-specifically describe the peakedness and/or heavy tail shape in a probability distribution via the standardized fourth moment [43]. For example, a Gaussian probability density function (pdf) has a kurtosis value of 3, whereas the hyperbolic secant pdf has a kurtosis value of 5 as its shape is both more peaked and heavier-tailed than the Gaussian shape. A convenient index called excess kurtosis is defined as the difference between the estimated kurtosis value for a given distribution with the value of the Gaussian pdf (i.e., 3), such that the excess kurtosis of a Gaussian pdf is 0. In applications of diffusion MRI, approaches have been developed to expand the logarithm of the exponential signal decay in the form of a Taylor series to arrive at higher order diffusion models [20, 44, 45]. Specifically, in Jensen et al. [20], the concept of excess kurtosis was connected to the Taylor series expansion of *b*, which, in its simplest scalar form, is given by,

$$S/S0=exp(-bD+\frac{1}{6}{b}^{2}{D}^{2}{K}_{\mathit{app}}),$$

(24)

where *K _{app}* is the apparent excess kurtosis. When

$${K}_{t}\equiv \frac{\langle {x}^{4}\rangle}{{\langle {x}^{2}\rangle}^{2}}-3.$$

(25)

Due to the parabolic form of the argument in Equation (24), there is a limit on the maximum *b*-value which can be fitted using this model. For typical diffusion MRI measurements in the human brain, the maximum values range from *b* ≤ 2000 – 4000 s/mm^{2} in order for the fitting function to monotonically decrease with increased diffusion weighting, as detailed in Jensen and Helpern [21]. In contrast to the parabolic form of Equations (24), (3), (8), (10), and (21) monotonically decrease as *b*→∞. However, Equations (3), (10), and (21) cannot be analytically expanded at *q* = 0 when *β* < 2 and, so, the second moment diverges for the diffusion propagator, *P*(*x*, *t*), which by extension, also applies to the stretched exponential models (about *q* or *b*) in Bennett et al. [14], Hall and Barrick [15], Magin et al. [16], Palombo et al. [17]. Although the second moment for these models diverge, it has been shown that rescaling methods lead to a pseudo MSD that is proportional to a trajectory that is faster than the Gaussian case of linear time dependence [25, 26]. On the other hand, Equations (8) and (23) can be expanded as analytic functions of *q* in order to compute the higher order moments and estimate the excess kurtosis as shown below.

By first utilizing the simple form of the Laplace-Fourier solution to time-fractional subdiffusion in Equation (7), the MSD can be computed by taking the second derivative of Equation (7) with respect to *q* in the limit of *q* → 0 and then performing a Laplace inversion,

$$\langle {x}^{2}(t)\rangle ={\mathcal{L}}^{-1}\underset{q\to 0}{lim}\{\frac{{d}^{2}p(q,s)}{d{q}^{2}}\},$$

(26)

which gives,

$$\langle {x}^{2}(t)\rangle ={\mathcal{L}}^{-1}\{2{D}_{\alpha ,2}{s}^{-(\alpha +1)}\},$$

(27)

and utilizing the common Laplace and time domain transform pair,

$$\frac{{t}^{\alpha}}{\mathrm{\Gamma}(\alpha +1)}={\mathcal{L}}^{-1}\{{s}^{-(\alpha +1)}\},$$

(28)

yields the form of the MSD,

$$\langle {x}^{2}(t)\rangle =\frac{2{D}_{\alpha ,2}}{\mathrm{\Gamma}(\alpha +1)}{t}^{\alpha},$$

(29)

as reported in Metzler and Klafter [25]. Extending the formalism of operating in the Laplace domain, the fourth moment of Equation (7) can be expressed as,

$$\langle {x}^{4}(t)\rangle ={\mathcal{L}}^{-1}\{24{({D}_{\alpha ,2})}^{2}{s}^{-(2\alpha +1)}\},$$

(30)

and inverting into the time domain using the transform in Equation (28) gives,

$$\langle {x}^{4}(t)\rangle =\frac{24{({D}_{\alpha ,2})}^{2}}{\mathrm{\Gamma}(2\alpha +1)}{t}^{2\alpha}.$$

(31)

Inserting Equations (29) and (31) into Equation (25), the excess kurtosis of the MLF for time-fractional subdiffusion is,

$${K}_{\mathit{MLF}}=6\frac{{\mathrm{\Gamma}}^{2}(\alpha +1)}{\mathrm{\Gamma}(2\alpha +1)}-3.$$

(32)

When *α* = 1, the MLF becomes the monoexponential CF and *K _{MLF}* = 0. For 0 <

Plot of Equation (32) for the kurtosis, *K*_{MLF}, computed in the Mittag-Leffler representation of subdiffusion vs. time-fractional derivative, *α*.

Interestingly, the form in Equation (32) coincides with the universal scaling law derived in Goychuk et al. [46], lim

$$\underset{t\to \infty}{lim}\frac{\langle x{(t)}^{2}\rangle}{{\langle x(t)\rangle}^{2}}=2\frac{{\mathrm{\Gamma}}^{2}(\alpha +1)}{\mathrm{\Gamma}(2\alpha +1)}-1,$$

(33)

which gives the relationship of the mean particle position to its MSD for subdiffusive diffusion dynamics described by *α*. These authors show through numerical simulations that the relationship in Equation (33) holds, within statistical error, under varying values of applied forces and temperatures, as is clearly illustrated in Figure 2 of Goychuk et al. [46]. Furthermore, a similar result to Equation (33) was derived in He et al. [47], showing this relationship is also a measure of ergodicity breaking behavior. Therefore, it is possible that Equations (32) and (33) reflect fundamental properties of biological tissue structure, which are minimally impacted by temperature and pressure, for diffusion MRI studies.

In order to illustrate the diffusion dynamics specified in Equations (3), (4), (8), and (10), one-dimensional random walk simulations were performed using the R software environment [48]. These numerical simulations provide a context to demonstrate how the order of the fractional derivatives in Equation (1) impact the governing statistics in a random walk’s distributions of waiting times and jump lengths. To simulate the sample path as governed by the CF in Equation (4), a walk with independent and identically distributed (iid) Gaussian jump lengths and iid Gaussian waiting time increments produces the dynamics described by Brownian motion in Figure 2. To simulate the sample path as governed by the CF in Equation (8), a walk of iid Gaussian step lengths, but with iid power-law waiting time increments, produces the dynamics described by time-fractional subdiffusion in Figure 3 with *α =* 0.75. To simulate the sample path as governed by the CF in Equation (10), a walk of iid power-law jump lengths and iid Gaussian waiting time increments produces the dynamics described by space-fractional superdiffusion in Figure 4 with *β =* 1.5. To simulate the sample path as governed by the CF in Equation (3), in a walk with step lengths and waiting times that are both independently iid power-law produces the dynamics described by time- and space-fractional diffusion in Figure 5 with *α =* 0.75 *β =* 1.5. In order to compute the iid power law behavior for the waiting times were selected from the Pareto distribution,

One dimensional random walk simulation of Brownian motion given by a Gaussian distribution of the waiting time parameter (*α* = 1) and a Gaussian distribution of the jump length parameter (*β* = 2).

One dimensional random walk simulation of time-fractional subdiffusion given by a power-law distribution of the waiting time parameter (*α* = 0.75) and a Gaussian distribution of the jump length parameter (*β* = 2).

One dimensional random walk simulation of space-fractional superdiffusion given by a Gaussian distribution of the waiting time parameter (*α* = 1) and a power-law distribution of the jump length parameter (*β* = 1.5).

One dimensional random walk simulation of space- and time-fractional anomalous diffusion given by a power-law distribution of the waiting time parameter (*α* = 0.75) and a power-law distribution of the jump length parameter (*β* = 1.5).

$$F(t)={(\frac{t}{c})}^{-{\scriptstyle \frac{1}{\alpha}}},$$

(34)

and the zero-mean corrected Pareto distribution for step lengths in the positive and negative directions,

$$F(x)={(\frac{x}{c})}^{-{\scriptstyle \frac{1}{\beta}}}-\frac{\beta}{\beta -1}{c}^{{\scriptstyle \frac{1}{\beta}}},$$

(35)

where *c* is a constant chosen to be 1 for these simulations to readily demonstrate the square root relationship between distance and time for diffusion, and *t* and *x* have arbitrary units with *t* ≥ 1 and |*x*| ≥ 1. The R codes were adapted from the examples fully described in chapter 5 of Meerschaert and Sikorskii [34].

One patient with chronic ischemic stroke was scanned on a 3 Tesla Siemens Trio MRI scanner with a 16 channel transmit/receive head coil (SiemensMedical Solutions, Erlangen, Germany). The imaging protocol was approved by Institutional Review Board at Northwestern University. Diffusion-weighted spin-echo echo planar imaging (SE-EPI) experiments were performed with the following pulse sequence parameters: echo time *TE* = 102 ms, repetition time *TR* = 6 s, Δ = 41.2 ms, *δ* = 30.6 ms, *b*-values = 0, 500, 1000, 3000, 4000 s/mm^{2}, 3 orthogonal diffusion weighted directions, number of averages *NA* = 6, inplane voxel resolution = 2 × 2 mm, voxel thickness = 4mm, 20 axial slices, scan time ~ 6 min. The raw diffusion weighted data were corrected for Rician noise by estimating the variance (*σ*^{2}) in the signal intensity of the ventricle at each b-value, such that
${S}_{rn}=\sqrt{{S}^{2}-2{\sigma}^{2}}$, which is a limited approach in assuming the spatial distribution of the noise is homogeneous in this multi-channel acquisition [49, 50]. The Rician noise-corrected diffusion weighted images were skull-stripped utilizing the Brain Extraction Tool of the FMRIB Software Library [51]. All skull-stripped and Rician noise-corrected diffusion weighted images were co-registered to the *b =* 0 image space using Statistical Parametric Mapping (SPM8, Wellcome Department of Cognitive Neurology, London, UK, http://www.fil.ion.ucl.ac.uk/spm). Using the Levenberg-Marquardt minimization algorithm in Matlab (Mathworks, Natick, MA), the averages of the 3 diffusion weighted direction data were fitted on a voxel-wise basis to Equation (24) and Equation (22) with the MLF algorithm in Podlubny [52], Gorenflo et al. [53]. Following estimations of *D* and *α*, the excess kurtosis, *K _{MLF}* was computed using the conversion provided in Equation (32). Additionally, data from

The Monte Carlo simulations in Figures 2–5 can be interpreted as the random motion of a single particle governed by the statistical conditions imposed by the order of the fractional derivatives in the generalized diffusion equation in Equation (1). Qualitatively, it can be observed that the magnitude of the distances in the sample paths are on the order of square-root of the time duration, a characteristic which is particularly evident in the Brownian motion case of Figure 2. For the case of time-fractional subdiffusion in Figure 3, the sample path clearly evolves slower than the square-root of time, whereas for space-fractional superdiffusion in Figure 4 the particle diverges faster than the Brownian case over time. Interestingly, for the time- and space-fractional diffusion in Figure 5 in which the simulation conditions have a ratio of 2*α*/*β =* 1, the sample path distance roughly evolves with the square root of time, however, the trajectory is clearly different from the Brownian motion simulated in Figure 2. Experimentally, the case for Brownian motion in Figure 2 is readily distinguishable from the time- and space-fractional diffusion scenario in Figure 5 by inspection of the CFs with respect to the monoexponential form.

In Figure 3, the apparent subdiffusive behavior of the particle appears as Brownian motion within certain time intervals, but is also punctuated by long periods in which the particle makes little or no random movement as a consequence of the heavy-tailed likelihood to wait longer than in the Gaussian case. In the counter condition shown in Figure 4, the apparent superdiffusive behavior of the particle appears as Brownian motion for most of the time, but is interrupted by short intervals when large steps are made to displace the particle to a new non-local position (i.e., pseudo-transport) as a consequence of the heavy-tailed likelihood to jump further than the Gaussian case. By combining these iid power-law conditions in both distance and time as shown in Figure 5, it can be observed that the particle is experiencing: moments of apparent Brownian motion, periods of extended rest, and apparent non-local transport within short periods of time.

In the diffusion MRI experiments, Table 1 shows the estimations of *α*, *K _{MLF}*,

The one dimensional, single particle random walk simulations with varying statistical conditions provide an illustrative means to conceptualize the differences between Gaussian and non-Gaussian diffusion. Of course, in MRI measurements of diffusion in biological tissue, the spatial resolution of a voxel is mesoscopic and therefore contains populations of proton spins diffusing about a biological microstructure medium that has varying properties of heterogeneity and anisotropy. When the diffusion MRI experiment is designed to sample multiple length and/or time scales (i.e., *b*-values) in order to probe the tissue microstructure, it is reasonable to consider that the bulk diffusion properties are non-Gaussian due to the heterogeneity of the tissue, reflected by the presence of neurons, glial cells, and microvessels. The CTRW framework is rooted in the statistical properties of the diffusion process in order to offer a physical description of the random motion, and therefore, in the biological context, subdiffusion, for example, can be envisaged as proton spin populations trapped by the restrictions of neurons and glia (Figure 3) whereas superdiffusion, can be envisaged as proton spin populations carried by microvesicular transport (Figure 4). Typically, in MRI experiments, the diffusion weightings are selected such that *b* > 500 s/mm^{2} in order to minimize the influence of bulk flow in the diffusion measurement [42]. Therefore, in such experiments, it is appropriate to consider a model which allows for the sensitization of sub- and Gaussian diffusion behaviors, as proposed in Equation (22).

For diffusion MRI, it should be highlighted here that the direct relationship between the Fourier transform of the diffusion propagator and the diffusion weighted signal holds in the limit of the short pulse approximation (i.e., *δ* Δ). Preclinical imaging spectrometers are able to produce very small *δ*, however most clinical MRI systems are not equipped to meet this relationship, though new advances in gradient hardware have become available in the most advanced systems to approach *δ* Δ [23]. Therefore, the estimates of kurtosis in this study, using either Equations (24) and (32), approach the true kurtosis of the pdf for the diffusion propagator in the limit that *δ* Δ. In current standard clinical implementation of the diffusion-weighted SE-EPI experiment, MRI systems are programmable to select a particular *b*-value, however the specific values of *δ* and Δ are dependent on the magnitude of *b* as well as the available maximum gradient amplitude, *G*. As a best practice, one would want to tune *δ* and Δ, such that different time and length scales are probed in order to characterize the properties of the tissue microstructure. However, in this study, to achieve a maximum *b*-value of 4000 s/mm^{2}, the scanner software set Δ = 41.2ms, *δ* = 30.6ms, which gives an effective diffusion time of 31ms. At this effective diffusion time, the sampled *q*-values are sensitive to length scales from ~ 2.8 – 7.9μm. Considering water is diffusing at a rate of ~ 0.6 – 1.0 s/mm^{2} in the WM and GM, the effective net water displacement (
$\text{using}~\sqrt{2Dt}$) is ~ 6 – 8μm. Therefore, depending on the microstructural heterogeneity and tortuosity in the biological tissue, the choices made to sample *q–*space can result in Gaussian diffusion (due to spatial averaging) for small values of *q* or exploit the non-Gaussian dynamics if *q* is large enough. As a result, higher order modeling of the diffusion attenuation signal (i.e., *α*, *K _{MLF}*, and

Utilizing the moment expansion of the time-fractional form of the MLF in Equations (29) and (31) provides an intimate link between kurtosis and subdiffusion through the Γ function and *α*. However, this link does not mean that Equations (22) and (24) are interchangeable, as each mathematical approach is a means to estimate the true kurtosis of the pdf for the diffusion propagator. In Table 1, if we compare the estimated values for *K _{MLF}* and

The typical signal fits for WM, GM, and CSF voxels can be visualized in Figure 7. The parabolic form of Equation (24) is apparent in A, C, and E with the estimations in the WM and GM data fitted to the decreasing side of the parabola and the CSF data fitted to both the decreasing and increasing sides of the parabola. Clearly, for the *b* = 3000 and 4000 s/mm^{2} CSF data, the noise floor has been reached and both Equations (24) and (22) fitted the higher *b*-value data points poorly, however it appears that Equation (24) is more sensitive to converging to spuriously high kurtosis estimates in the ventricles and necrosed tissue due residual signal noise. Nevertheless, over a limited range of *b*-values, Equations (22) and (24) can produce similar, but not identical, estimations of kurtosis. As Equation (22) is not bounded by a maximum *b*-value, the MLF provides the opportunity to more completely sample *q*-space in order to more accurately estimate the true kurtosis of the diffusion propagator, which is an advantage as very high strength (300 mT/m) gradients have recently become available for diffusion imaging [23]. For example, by removing the *b* = 4000 s/mm^{2} data and fitting to Equation (24), estimates for *K _{app}* in the WMROI increase by 11.8±4.9% and, in the GM ROI,

We have presented new, simplified fitting forms for the MLF as a three-parameter model in Equation (21) (for potential subdiffusion and superdiffusion) and a two-parameter model in Equation (22) (for potential subdiffusion only). The concepts of subdiffusion, superdiffusion, and Brownian motion have been simulated to illustrate the physical consequences of the movement of a particle in the statistical context of the CTRW theory, which potentially can have biological correlates in diffusion MRI. We have computed the kurtosis (*K _{MLF}*) for time-fractional form of the MLF, which provides a context to relate subdiffusion (

**Funding**

This work has been funded by a grant from the Whitaker International Program of the Institute of International Education and by the National Institutes of Health grant NIDCD 1P50DC012283-01A1.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

1. Tanner JE, Stejskal EO. Restricted self-diffusion of protons in colloidal systems by the pulsed-gradient, spin-echo method. J Chem Phys. 1968;49:1768–77. doi: 10.1063/1.1670306. [Cross Ref]

2. Le Bihan D, Breton E, Lallemand D, Grenier P, Cabanis E, Laval-Jeantet M. MR imaging of intravoxel incoherent motions: application to diffusion and perfusion in neurologic disorders. Radiology. 1986;161:401–7. doi: 10.1148/radiology.161.2.3763909. [PubMed] [Cross Ref]

3. Holz M, Heil SR, Sacco A. Temperature-dependent self-diffusion coefficients of water and six selected molecular liquids for calibration in accurate 1HNMR PFG measurements. 2000;2:4740–2. doi: 10.1039/B005319H. [Cross Ref]

4. Norris DG, Niendorf T, Leibfritz D. Healthy and infarcted brain tissues studied at short diffusion times: the origins of apparent restriction and the reduction in apparent diffusion coefficient. NMR Biomed. 1994;7:304–10. doi: 10.1002/nbm.1940070703. [PubMed] [Cross Ref]

5. van Gelderen P, de Vleeschouwer MH, DesPres D, Pekar J, van Zijl PC, Moonen CT. Water diffusion and acute stroke. Magn Reson Med. 1994;31:154–63. doi: 10.1002/mrm.1910310209. [PubMed] [Cross Ref]

6. Niendorf T, Dijkhuizen RM, Norris DG, van Lookeren Campagne M, Nicolay K. Biexponential diffusion attenuation in various states of brain tissue: implications for diffusion-weighted imaging. Magn Reson Med. 1996;36:847–57. doi: 10.1002/mrm.1910360607. [PubMed] [Cross Ref]

7. Inglis BA, Bossart EL, Buckley DL, Wirth ED, Mareci TH. Visualization of neural tissue water compartments using biexponential diffusion tensor MRI. Magn Reson Med. 2001;45:580–87. doi: 10.1002/mrm.1079. [PubMed] [Cross Ref]

8. Thelwall PE, Grant SC, Stanisz GJ, Blackband SJ. Human erythrocyte ghosts: exploring the origins of multiexponential water diffusion in a model biological tissue with magnetic resonance. Magn Reson Med. 2002;48:649–57. doi: 10.1002/mrm.10270. [PubMed] [Cross Ref]

9. Assaf Y, Basser PJ. Composite hindered and restricted model of diffusion (CHARMED)MR imaging of the human brain. Neuroimage. 2005;27:48–58. doi: 10.1016/j.neuroimage.2005.03.042. [PubMed] [Cross Ref]

10. Assaf Y, Blumenfeld-Katzir T, Yovel Y, Basser PJ. AxCaliber: a method for measuring axon diameter distribution from diffusion MRI. Magn Reson Med. 2008;59:1347–54. doi: 10.1002/mrm.21577. [PMC free article] [PubMed] [Cross Ref]

11. Alexander DC, Hubbard PL, Hall MG, Moore EA, Ptito M, Parker GJM, et al. Orientationally invariant indices of axon diameter and density from diffusion MRI. Neuroimage. 2010;52:1374–89. doi: 10.1016/j.neuroimage.2010.05.043. [PubMed] [Cross Ref]

12. Zhang H, Schneider T, Wheeler-Kingshott CA, Alexander DC. NODDI: practical *in vivo* neurite orientation dispersion and density imaging of the human brain. Neuroimage. 2012;61:1000–16. doi: 10.1016/j.neuroimage.2012.03.072. [PubMed] [Cross Ref]

13. Panagiotaki E, Walker-Samuel S, Siow B, Johnson SP, Rajkumar V, Pedley RB, et al. Noninvasive quantification of solid tumor microstructure using VERDICT MRI. Cancer Res. 2014;74:1902–12. doi: 10.1158/0008-5472.CAN-13-2511. [PubMed] [Cross Ref]

14. Bennett KM, Schmainda KM, Bennett RT, Rowe DB, Lu H, Hyde JS. Characterization of continuously distributed cortical water diffusion rates with a stretched-exponential model. Magn Reson Med. 2003;50:727–34. doi: 10.1002/mrm.10581. [PubMed] [Cross Ref]

15. Hall MG, Barrick TR. From diffusion-weighted MRI to anomalous diffusion imaging. Magn Reson Med. 2008;59:447–55. doi: 10.1002/mrm.21453. [PubMed] [Cross Ref]

16. Magin RL, Abdullah O, Baleanu D, Zhou XJ. Anomalous diffusion expressed through fractional order differential operators in the Bloch-Torrey equation. J Magn Reson. 2008;190:255–70. doi: 10.1016/j.jmr.2007.11.007. [PubMed] [Cross Ref]

17. Palombo M, Gabrielli A, De Santis S, Cametti C, Ruocco G, Capuani S. Spatiotemporal anomalous diffusion in heterogeneous media by nuclear magnetic resonance. J Chem Phys. 2011;135:34504. doi: 10.1063/1.3610367. [PubMed] [Cross Ref]

18. Özarslan E, Basser PJ, Shepherd TM, Thelwall PE, Vemuri BC, Blackband SJ. Characterization of anomalous diffusion from MR signal may be a new probe to tissue microstructure. Conf Proc IEEE Eng Med Biol Soc. 2006;1:2256–9. doi: 10.1109/IEMBS.2006.259651. [PubMed] [Cross Ref]

19. Özarslan E, Shepherd TM, Koay CG, Blackband SJ, Basser PJ. Temporal scaling characteristics of diffusion as a new MRI contrast: findings in rat hippocampus. Neuroimage. 2012;60:1380–93. doi: 10.1016/j.neuroimage.2012.01.105. [PMC free article] [PubMed] [Cross Ref]

20. Jensen JH, Helpern JA, Ramani A, Lu H, Kaczynski K. Diffusional kurtosis imaging: the quantification of non-gaussian water diffusion by means of magnetic resonance imaging. Magn Reson Med. 2005;53:1432–40. doi: 10.1002/mrm.20508. [PubMed] [Cross Ref]

21. Jensen JH, Helpern JA. MRI quantification of non-Gaussian water diffusion by kurtosis analysis. NMR Biomed. 2010;23:698–710. doi: 10.1002/nbm.1518. [PMC free article] [PubMed] [Cross Ref]

22. Wu EX, Cheung MM. MR diffusion kurtosis imaging for neural tissue characterization. 2010;23:836–48. doi: 10.1002/nbm.1506. [PubMed] [Cross Ref]

23. Van Essen DC, Ugurbil K, Auerbach E, Barch D, Behrens TEJ, Bucholz R, et al. The human connectome project: a data acquisition perspective. Neuroimage. 2012;62:2222–31. doi: 10.1016/j.neuroimage.2012.02.018. [PMC free article] [PubMed] [Cross Ref]

24. Ingo C, Magin RL, Colon-Perez L, Triplett W, Mareci TH. On random walks and entropy in diffusion-weighted magnetic resonance imaging studies of neural tissue. Magn Reson Med. 2014;71:617–27. doi: 10.1002/mrm.24706. [PMC free article] [PubMed] [Cross Ref]

25. Metzler R, Klafter J. The random walk’s guide to anomalous diffusion: a fractional dynamics approach. Phys Rep. 2000;339:1–77. doi: 10.1016/S0370-1573(00)00070-3. [Cross Ref]

26. Metzler R, Jeon JH, Cherstvy AG, Barkai E. Anomalous diffusion models and their properties: non-stationarity, non-ergodicity, and ageing at the centenary of single particle tracking. Phys Chem Chem Phys. 2014;16:24128–64. doi: 10.1039/C4CP03465A. [PubMed] [Cross Ref]

27. Capuani S, Palombo M, Gabrielli A, Orlandi A, Maraviglia B, Pastore FS. Spatio-temporal anomalous diffusion imaging: results in controlled phantoms and in excised human meningiomas. Magn Reson Imaging. 2013;31:359–65. doi: 10.1016/j.mri.2012.08.012. [PubMed] [Cross Ref]

28. Gadelkarim JJ, Magin RL, Meerschaert MM, Capuani S, Palombo M, Kumar A, et al. Fractional order generalization of anomalous diffusion as a multidimensional extension of the transmission line equation. IEEE J Emerging Sel Top Circ Syst. 2013;3:432–41. doi: 10.1109/JETCAS.2013.2265795. [Cross Ref]

29. Gorenflo R, Vivoli A, Mainardi F. Discrete and continuous random walk models for space-time fractional diffusion. Nonlin Dyn. 2004;38:101–16. doi: 10.1007/s11071-004-3749-5. [Cross Ref]

30. Brockmann D, Hufnagel L, Geisel T. The scaling laws of human travel. Nature. 2006;439:462–5. doi: 10.1038/nature04292. [PubMed] [Cross Ref]

31. Metzler R, Glackle WG, Nonnenmacher TF. Fractional model equation for anomalous diffusion. Physica A Stat Mech Appl. 1994;211:13–24. doi: 10.1016/0378-4371(94)90064-7. [Cross Ref]

32. Magin RL. Fractional Calculus in Bioengineering. Connecticut: Begell House Publishers; 2006.

33. Zaslavsky GM. Hamiltonian Chaos and Fractional Dynamics. New York, NY: Oxford University Press; 2005.

34. Meerschaert MM, Sikorskii A. Stochastic Models for Fractional Calculus. Vol. 43. Berlin: De Gruyter; 2012.

35. Podlubny I. Fractional Differential Equations. San Diego: Academic Press; 1999.

36. Mittag-Leffler GM. Sur la nouvelle fonction E*α* (x) CR Acad Sci Paris. 1903;137:554–8.

37. Mittag-Leffler G. Sur la representation analytique d’une branche uniforme d’une fonction monogene. Acta Math. 1905;29:101–81. doi: 10.1007/BF02403200. [Cross Ref]

38. Haubold HJ, Mathai AM, Saxena RK. Mittag-Leffler functions and their applications. J Appl Math. 2011;2011:298628.

39. Wheatcraft SW, Meerschaert MM. Fractional conservation of mass. Adv Water Resour. 2008;31:1377–81. doi: 10.1016/j.advwatres.2008.07.004. [Cross Ref]

40. Leszczynski JS. An Introduction to Fractional Mechanics. Częstochowa: Częstochowa University of Technology; 2011.

41. West BJ, Bologna M, Grigolini P. Physics of Fractal Operators. New York, NY: Springer; 2003.

42. Le Bihan D, Turner R, Douek P, Patronas N. Diffusion MR imaging: clinical applications. Am J Roentgenol. 1992;159:591–9. doi: 10.2214/ajr.159.3.1503032. [PubMed] [Cross Ref]

43. Balanda KP, MacGillivray HL. Kurtosis: a critical review. Am Stat. 1988;42:111–9.

44. Liu C, Bammer R, Acar B, Moseley ME. Characterizing non-gaussian diffusion by using generalized diffusion tensors. Magn Reson Med. 2004;51:924–37. doi: 10.1002/mrm.20071. [PubMed] [Cross Ref]

45. Yablonskiy DA, Bretthorst GL, Ackerman JJH. Statistical model for diffusion attenuated MR signal. Magn Reson Med. 2003;50:664–9. doi: 10.1002/mrm.10578. [PMC free article] [PubMed] [Cross Ref]

46. Goychuk I, Heinsalu E, Patriarca M, Schmid G, Hänggi P. Current and universal scaling in anomalous transport. Phys Rev E. 2006;73:020101. doi: 10.1103/PhysRevE.73.020101. [PubMed] [Cross Ref]

47. He Y, Burov S, Metzler R, Barkai E. Random time-scale invariant diffusion and transport coefficients. Phys Rev Lett. 2008;101:058101. [PubMed]

48. Team RC. Technical Report. R Foundation for Statistical Computing; Vienna: 2012. R: A Language and Environment for Statistical Computing.

49. McGibney G, Smith MR. An unbiased signal-to-noise ratio measure for magnetic resonance images. Med Phys. 1993;20:1077–78. doi: 10.1118/1.597004. [PubMed] [Cross Ref]

50. Miller AJ, Joseph PM. The use of power images to perform quantitative analysis on low SNR MR images. 1993;11:1051–6. [PubMed]

51. Smith SM. Fast robust automated brain extraction. Hum Brain Mapp. 2002;17:143–55. doi: 10.1002/hbm.10062. [PubMed] [Cross Ref]

52. Podlubny I. The Mittag-Leffler Function. 2009 Available online at: www.mathworks.com/matlabcentral/fileexchange/8738.

53. Gorenflo R, Loutchko J, Luchko Y. E*α*, *β* (z) and its derivative. Fract Calculus Appl Anal. 2002;5:491–518.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |