Anisotropic fractional diffusion tensor imaging

1. Introduction

The structural complexity of the human brain is manifest at each level of functional organization: synapses, axons, neurons, cortical layers, fiber tracts, and cerebral convolutions (gyri and sulci) (Schaltenbrand and Wahren, 1998). Magnetic resonance imaging (MRI) in general, and diffusion tensor imaging (DTI) in particular, exhibit contrast that reflects tissue heterogeneity and anisotropy in both the white and the gray matter (Mori, 2006). The overall goal of these imaging modalities is to provide spatial maps of structural features that correspond to the specific neural networks that provide the basis for sensory awareness, memory, cognition and coordinated movement (Le Bihan, 1995). Disruption of these neural pathways is a hallmark of trauma, stroke, tumors and degenerative disease. Although MRI and DTI are useful clinical tools for diagnosis and treatment monitoring, their typical voxel resolution (1 mm³) is orders of magnitude above that of a single cell (10 µm³) (Johansen-Berg and Behrens, 2009). Therefore there is a need to probe sub-voxel structure to improve both the sensitivity and the specificity of diagnosis. Since water movement within the voxel leads to MR signal attenuation that reflects collisions with molecules, membranes, and axonal fibers (Haacke et al., 1999) we anticipate that stochastic models of diffusion (isotropic, anisotropic, restricted, hindered, Gaussian, nonGaussian) can be used to encode sub-millimeter structure.

2. Fractional DTI

The connection between diffusion and magnetic resonance for water protons is described by the Bloch–Torrey equation (Torrey, 1956; Haacke et al., 1999; Callaghan, 2011). Solving the Bloch–Torrey equation for an anisotropic material, such as brain white matter (WM), provides the basis for DTI (Le Bihan, 1995). In standard DTI, a pair of trapezoidal gradient pulses is added to the MR imaging sequence (Mori, 2006). The acquired diffusion-weighted (DW) signal S decays in a manner dependent upon the diffusion gradient strength, G, gradient duration, δ, and the time interval, Δ, between gradient pulses. The resultant decay can be modeled (Haacke et al., 1999) by the equation

S = S 0 exp ( - ( γ G δ ) 2 g · D g · ( Δ - δ / 3 ) )

(2.1)

S ( b , g ) = S 0 exp ( - b g · D g )

(2.2)

∂ t p ( x , t ) = ∇ · D ∇ p ( x , t )

(2.3)

Since D is symmetric and positive definite, there is an orthonormal basis of eigenvectors v ₁, … ,  v _d with corresponding eigenvalues a_i such that D v _i  = a_i v _i for 1 ≤ i ≤ d. For any k  ∈ ℝ ^d we can write k  =  ∑ j = 1 d k j v j where k_j =  k  ·  v _j . Then v _i  · D v _j  = 0 if j ≠ i and v _i  · D v _i  = a_i. It follows easily that

p ∧ ( k , t ) = exp [ - t ∑ i = 1 d a i k i 2 ]

(2.4)

D = ( 1 / 2 0 0 2 )

(2.5)

OPEN IN VIEWER

Figure 1.

Level sets of the Fourier solution (2.4) to the traditional diffusion equation in d = 2 dimensions with diffusion tensor (2.5) at b = 1.

In many applications (Metzler and Klafter, 2000, 2004; Mainardi, 2010; Meerschaert and Sikorskii, 2012), a diffusing front spreads at a different rate than the t^{1 / 2} predicted by the traditional diffusion equation. This can be captured by introducing fractional derivatives into the diffusion model. For simplicity, let us focus on the isotropic diffusion model where D = DI for some positive constant D, and where I is the d × d identity matrix. Then the diffusion equation (2.3) reduces to ∂_tp( x , t) = DΔp( x , t), and its point source solution has Fourier transform ( k , t) = exp(−Dt‖ k ‖²) where ‖ k ‖² =  k  ·  k . The level sets of the solution p( x , t) are circles in two dimensions, or spheres in three dimensions.

The fractional Laplacian is an isotropic space-fractional derivative, defined so that Δ ^{α/ 2}f( x ) has Fourier transform −‖ k ‖ ^α f ∧ ( k ) with 0 < α < 2. This reduces to the traditional Laplacian when α = 2. The isotropic space-fractional diffusion equation

∂ t p ( x , t ) = D Δ α / 2 p ( x , t )

(2.6)

3. Fractional Bloch–Torrey equation

The traditional Bloch–Torrey equation

∂ t S = - ( i γ x · G ) S + ∇ · D ∇ S

(3.1)

L : = γ ∫ 0 t G ( τ ) d τ

A ' A S - ( i x · γ G ) S = - i γ ( x · G ) S + ∇ · D ∇ S

A ( t ) = exp [ - ∫ 0 t L ( τ ) · D L ( τ ) d τ ]

The simplest space-fractional Bloch–Torrey equation

∂ t S = - ( i γ x · G ) S + D 0 Δ α / 2 S

(3.2)

A ' A S = D 0 Δ α / 2 S

A ( t ) = exp [ - D 0 ∫ 0 t ‖ L ( τ ) ‖ α d τ ]

S = S 0 exp [ - D 0 ( γ G δ ) α ( Δ - α - 1 α + 1 δ ) ]

(3.3)

S = S 0 exp ( - b α D )

(3.4)

4. The need for anisotropic fractional DTI models

In MRI experiments, it is often observed (Bennett et al., 2003; Hall and Barrick, 2008; Ingo et al., 2014) that the acquired DW signal S follows the stretched exponential model (3.4) for high b values. In applications to brain imaging, it is also found that the parameter α varies with direction: an indication of anisotropy (Hall and Barrick, 2012). For example, in one experiment (Ingo et al., 2014), formalin-fixed brains from normal, adult rats were soaked in Fluorinert to reduce magnetic susceptibility and imaged ex vivo in a Bruker 750 MHz spectrometer (17.6 T, 89 mm bore). A pulsed gradient stimulated echo diffusion sequence was used with pulse repetition time of 2 s, echo time of 28 ms, in-plane resolution of 190 µm and slice thickness of 1 mm. The signal S was acquired in six different vector directions g ₁ = (0, 0, 1) ^T , g ₂ = (0.89, 0, 0.45) ^T , g ₃ = (0.28, 0.09, 0.45) ^T , g ₄ = (−0.72, −0.53, 0.45) ^T , g ₅ = (0.28, −0.85, 0.45) ^T , and g ₆ = (−0.72, −0.53, 0.45) ^T with 10 different b values ranging up to a maximum value of 26,190 s/mm². In this experiment, Δ (17.5 ms) and δ (3.5 ms) were kept constant and G was scaled to increase with the b value. Under these conditions, the short-pulse condition δ ≪ Δ holds.

Next, we validate the stretched exponential model (3.4) using linear regression. Taking logs in the model yields ln(S/S₀) = −Db^α, and taking logs again produces

s : = ln ( - ln ( S / S 0 ) ) = ln D + α ln b

(4.1)

OPEN IN VIEWER

Figure 2.

Plot of s versus ln b in six different directions, to validate the stretched exponential model (4.1).

To illustrate the anisotropic nature of DTI, we compare the slopes α from the straight line fit of s versus ln b for each direction j = 1, 2, … , 6. The results are summarized in Table 1. It is clear that the power law slope depends on direction. For example, the slope for direction 4 (corresponding to direction vector g ₄) is α = 0.379 ±0.023 which is significantly different from the α = 0.294 value in direction 2. To obtain a formal confidence interval for these α values, one can use the standard t-interval from linear regression theory. Since the sample size is n = 7, the 95% confidence interval is α ± 2.571 SE where α is the estimate in the first row of Table 1, SE is the standard error in the second row of Table 1, and 2.571 is the 97.5th percentile of the standard t distribution with n − 2 = 5 degrees of freedom. For example, we are 95% confident that the correct α value for direction 1 lies in the interval (0.300, 0.336).

OPEN IN VIEWER

Table 1.

Best-fit α values via linear regression on the data in Figure 2 for six different directions, demonstrating anisotropy.

Direction	1	2	3	4	5	6
α	0.318	0.294	0.303	0.379	0.305	0.362
Standard error	0.007	0.014	0.022	0.023	0.018	0.018

Previous work has shown α to be a biomarker that reflects tissue heterogeneity (Bennett et al., 2003; Ozarslan and Mareci, 2003; Hall and Barrick, 2008; Magin et al., 2008; Zhou et al., 2010; Palombo et al., 2011; Ingo et al., 2014; Magin et al., 2014). In particular, the stretched exponential parameter α exhibits a lower value in more tortuous porous materials, and more heterogeneous tissue. Since the heterogeneity parameter α also varies with direction, it would be advantageous to include anisotropy in the fractional DTI model, to capture this effect.

5. Anisotropic fractional diffusion models for DTI

In the previous section, we demonstrated that the anisotropy parameter α in a DTI model can vary with direction. An open challenge in MRI theory is to develop a suitable anisotropic model for fractional diffusion that captures this anisotropy, and can be fit to data acquired within the constraints of a clinical MR scan. The scan generates multidimensional data that must be quickly and accurately presented to the radiologist for analysis. Fractional order models are attractive because they capture tissue complexity in a small set of parameters. Within the continuous time random walk (CTRW) paradigm, for example, the random motion of water is hindered or restricted by tissue heterogeneity that alters the waiting times and jump increments in a manner simply expressed by power laws (Metzler and Klafter, 2000, 2004; Meerschaert and Sikorskii, 2012). In the remainder of this paper, we survey existing anisotropic models for fractional DTI, and also discuss some potential new models.

5.1. Anisotropic fractional diffusion

A recent paper of Hanyga and Magin (2014) proposed a new space-fractional diffusion model that seems well-suited for applications to DTI. The model is

∂ t p ( x , t ) = Q p ( x , t )

(5.1)

where the anisotropic fractional derivative operator Q is defined in terms of Fourier transforms. Define

Q ∧ ( k ) = - ∫ y ∈ Ω | y · k | α ( y ) m ( d y )

(5.2)

where Ω = {y ∈ ℝ ^d : ‖y‖ = 1} is the unit sphere in d-dimensional Euclidean space, m(dy) is a finite Borel measure on the sphere, and α(y) is a symmetric function α(y) = α(−y) on the sphere that takes values in the interval (0, 2). Then one defines Qf(x) to be the function with Fourier transform Q ∧ (k)f^( k ). Hanyga and Magin (2014) continue to prove that solutions to (5.1) remain nonnegative for a nonnegative initial condition p( x , 0) = f( x ) ≥ 0. Next we provide an alternative proof of this fact, by showing that the solutions to (5.1) are the probability densities of a Lévy process (Sato, 1999; Meerschaert and Sikorskii, 2012, Section 4.3).

Proposition 5.1

There exists a Lévy process X(t) that satisfies

p ∧ ( k , t ) = E [ e - i k · X ( t ) ] = e t Q ∧ ( k )

(5.3)

for all k ∈ ℝ ^d and all t ≥ 0, for any symmetric index function α: Ω → (0, 2) and any finite Borel measure m(dy) on the unit sphere.

Proof

Any Lévy process {X(t) : t ≥ 0} is determined by the distribution of X: = X(1), which can be specified using the Lévy representation (Meerschaert and Scheffler, 2001, Theorem 3.1.11): a random vector X on ℝ ^d is infinitely divisible if and only if we can write E[e^{− i k  ·  X} ] = exp( Q ∧ (k)), where

Q ∧ ( k ) = - i a · k + 12 k · A k - ∫ x ≠ 0 ( e i k · x - 1 - i k · x 1 + ‖ x ‖ 2 ) φ ( d x )

(5.4)

for a ∈ ℝ ^d , A a nonnegative definite d × d matrix, and φ a σ-finite Borel measure on ℝ ^d \{0} such that

∫ x ≠ 0 min { 1 , ‖ x ‖ 2 } φ ( d x ) < ∞

(5.5)

The triple [a, A, φ] is unique. Next we note that, in the one-dimensional case d = 1, there exists an infinitely divisible random variable X such that E[e^{− ikX} ] = exp(−|k| ^α ) for any 0 < α < 2, and in this case, it follows from Meerschaert and Scheffler (2001, Lemma 7.3.10) (for 0 < α < 1), Meerschaert and Scheffler (2001, Lemma 7.3.11) (for 1 < α < 2), and Meerschaert and Scheffler (2001, Lemma 7.3.12) (for α = 1) that this random variable has Lévy representation [0, 0, φ_α] where

φ α ( d x ) = C α 2 α | x | - α - 1 d x

(5.6)

with

C α = 1 - α Γ ( 2 - α ) cos ( π α / 2 ) for 0 < α < 1 or 1 < α < 2

(5.7)

and C₁ = 2/π. Then we have

| k | α = ∫ x ≠ 0 ( e ikx - 1 - ikx 1 + x 2 ) φ α ( d x )

(5.8)

for each 0 < α < 2. Next, define an infinitely divisible random vector X on ℝ ^d (e.g. let d = 3) by specifying its Lévy representation [0, 0, φ] where x = ry in polar coordinates r > 0 and ‖y‖ = 1, and

φ ( d x ) = φ ( dr , d y ) = C α ( y ) α ( y ) r - α ( y ) - 1 dr m ¯ ( d y )

(5.9)

where m(dy) = [m(dy) + m(−dy)]/2 is the symmetrized version of the measure m(dy). Then the random vector X has Lévy representation E[e^{− i k  ·  X} ] = exp( Q ∧ (k)), where

- Q ∧ ( k ) = ∫ y ∈ Ω ∫ 0 ∞ ( e ir k · y - 1 - ir k · y 1 + r 2 ( k · y ) 2 ) × C α ( y ) α ( y ) r - α ( y ) - 1 dr m ¯ ( d y ) = ∫ y ∈ Ω ∫ 0 ∞ ( e ir k · y - 1 - ir k · y 1 + r 2 ( k · y ) 2 ) × C α ( y ) 2 α ( y ) r - α ( y ) - 1 dr m ( d y ) + ∫ y ∈ Ω ∫ 0 ∞ ( e ir k · y - 1 - ir k · y 1 + r 2 ( k · y ) 2 ) × C α ( y ) 2 α ( y ) r - α ( y ) - 1 dr m ( - d y ) = ∫ y ∈ Ω ∫ 0 ∞ ( e ir k · y - 1 - ir k · y 1 + r 2 ( k · y ) 2 ) × C α ( y ) 2 α ( y ) r - α ( y ) - 1 dr m ( d y ) + ∫ y ∈ Ω ∫ 0 ∞ ( e - ir k · y - 1 - - ir k · y 1 + r 2 ( k · y ) 2 ) × C α ( y ) 2 α ( y ) r - α ( y ) - 1 dr m ( d y ) = ∫ y ∈ Ω ∫ 0 ∞ ( e ir k · y - 1 - ir k · y 1 + r 2 ( k · y ) 2 ) × C α ( y ) 2 α ( y ) r - α ( y ) - 1 dr m ( d y ) + ∫ y ∈ Ω ∫ - ∞ 0 ( e ir k · y - 1 - ir k · y 1 + r 2 ( k · y ) 2 ) × C α ( y ) 2 α ( y ) | r | - α ( y ) - 1 dr m ( d y ) = ∫ y ∈ Ω ∫ r ≠ 0 ( e ir ( k · y ) - 1 - ir ( k · y ) 1 + r 2 ) × C α ( y ) 2 α ( y ) | r | - α ( y ) - 1 dr m ( d y ) + ∫ y ∈ Ω ∫ r ≠ 0 ( ir ( k · y ) 1 + r 2 - ir ( k · y ) 1 + r 2 ( k · y ) 2 ) × C α ( y ) 2 α ( y ) | r | - α ( y ) - 1 dr m ( d y ) = ∫ y ∈ Ω | k · y | α ( y ) m ( d y )

(5.10)

in view of (5.8), since the integral in the next to last line equals zero by symmetry. This shows that (5.3) holds. Since C_α is a bounded continuous function on the interval 0 < α < 2, it follows easily that (5.5) holds, so (5.9) is a Lévy measure.□

We say that a random vector X is full if it is not supported on a lower-dimensional hyperplane, that is, if there is no unit vector y ∈ Ω such that X · y = 0 with probability one.

Proposition 5.2

If ∫ _{y   ∈   Ω}|y · w|^{α( y )}m(dy) > 0 for every w ∈ Ω, and α(y) ≥ α₀ > 0 for all y ∈ Ω, then the infinitely divisible random variable X(t) in Proposition 5.1 is full, and has a density p(x, t) with respect to Lebesgue measure for any t > 0.

Proof

Define p ∧ ( k , t ) = e t Q ∧ ( k ) for all k ∈ ℝ ^d and all t ≥ 0. The Fourier inversion theorem (Meerschaert and Scheffler, 2001, Theorem 1.3.7) implies that

p ( x , t ) = ( 2 π ) - d ∫ e i k · x p ∧ ( k , t ) d k

(5.11)

is the function with Fourier transform p ∧ (k, t), so long as the integral ∫| p ∧ (k, t)|dk < ∞. Since Q ∧ (k) ≤ 0 for any k ∈ ℝ ^d , it follows that 0 ≤  p ∧ (k, t) ≤ 1 for all k ∈ ℝ ^d and all t ≥ 0. Hence it suffices to check that ∫_{‖ k ‖ ≥ 1}| p ∧ (k, t)|dk < ∞. Adopt the polar coordinates k = ρw where ρ > 0 and ‖w‖ = 1. Since α(y) ≥ α₀ > 0 for all y ∈ Ω, we have

- Q ∧ ( k ) = ∫ y ∈ Ω ρ α ( y ) | y · w | α ( y ) m ( d y ) ≥ ρ α 0 ∫ y ∈ Ω | y · w | α ( y ) m ( d y )

for all ‖k‖ ≥ 1. It follows from the Dominated Convergence Theorem (Rudin, 1976, Theorem 11.32) that g(w) := ∫ _{y ∈Ω}|y · w| ^{α ( y )}m(dy) is a continuous function on the compact set Ω, and since g(w) > 0 by assumption, it follows that g(w) ≥ g₀ > 0 for all w ∈ Ω. Then − Q ̂ ( k ) ≥ g 0 ρ α 0 for all ‖k‖ ≥ 1. Then we have

∫ ‖ k ‖ ≥ 1 | p ∧ ( k , t ) | d k ≤ ∫ w ∈ Ω ∫ 0 ∞ e - g 0 ρ α 0 d ρ m ( d w ) ≤ C 0 m ( Ω )

where C 0 = ∫ 0 ∞ e - g 0 ρ α 0 d ρ < ∞ . This shows that (5.11) holds. If X(t) were not full, then we would have X(t) · w = 0 for some unit vector w, and then we would have E(e^{− i w  ·  X ( t )}) = e^{t Q ∧ ( w )} = 1, hence Q ∧ (w) = 0. But this contradicts − Q ∧ (k) ≥ g₀ρ^α₀, and so X(t) is full for every t > 0.□

Remark 5.3

Since p(x, t) is a probability density for any t > 0, it follows that p(x, t) ≥ 0 for all x ∈ ℝ ^d and all t > 0. With some additional work, it should be possible to show that p(x, t) > 0 for all x ∈ ℝ ^d and all t > 0.

Remark 5.4

It should not be hard to extend these arguments to the asymmetric case

Q ∧ ( k ) = ∫ y ∈ Ω ( y · i k ) α ( y ) m ( d y )

(5.12)

which reduces to the case (5.2) if the measure m(dy) is symmetric, that is, m(dy) = m(−dy). One just has to be a bit careful about the centering constants.

Remark 5.5

A variety of stable-like processes have been considered in the literature, but the process constructed in Proposition 5.1 seems new. Bass (1988) considers a stable-like process in one dimension with jump intensity

φ ( x , dy ) = C α ( x ) 2 α | y | - α ( x ) - 1 dy

which behaves locally like an α-stable process whose index varies in space. Bass and Tang (2009) consider a d-dimensional stable-like process with jump intensity φ(x, dy) = A(x,y)‖y‖^{− α −d}dy where A(x,y) > 0 is bounded away from zero and infinity. That model exhibits mild anisotropy, as opposed to the strong anisotropy in the Hanyga model. It would certainly be interesting to explore the mathematical properties of the Lévy process in Proposition 5.1 in more detail.

5.2. Anisotropic fractional Bloch–Torrey equation

Here we propose a new anisotropic fractional Bloch–Torrey equation

∂ t S = - ( i γ x · G ) S + D 0 Q S

(5.13)

which can be solved by the method introduced in Section 3. Assume the same solution form S = S₀Ae^{− i x  ·  L} as before, and compute

A ' A S = D 0 Q

Next compute Q ∧ S using Fourier transforms. Recall that the function f( x ) = e^{− i a  ·  x} has Fourier transform f ∧ ( k )  =(2π)^dδ( k  +  a ). Then note that D₀ Q ∧ S has Fourier transform

D 0 Q ∧ ( k ) S ∧ ( k , t ) = - D 0 ∫ y ∈ Ω | y · k | α ( y ) m ( d y ) S ∧ ( k , t )

Inverting as in Section 3, it follows that

A ( t ) = exp [ - D 0 ∫ 0 t ∫ y ∈ Ω | y · L ( τ ) | α ( y ) m ( d y ) d τ ]

in this case. For a Stejskal–Tanner pulse sequence, the solution reduces to

S = S 0 exp [ - D 0 ∫ y ∈ Ω | γ G δ · y | α ( y ) ( Δ - α ( y ) - 1 α ( y ) + 1 δ ) m ( d y )

(5.14)

If the mixing measure m(d y ) is concentrated on d point masses on an arbitrary set of coordinate axes v ₁, … ,  v _d which need not be orthogonal, this reduces to a model recently proposed and tested by GadElkarim et al. (2013). That model has the solution

S = S 0 exp [ - ∑ j = 1 d D j | γ G δ · v j | α j ( Δ - α j - 1 α j + 1 δ ) ]

(5.15)

where D_j = D₀m( v _j ), which agrees with GadElkarim et al. (2013, equation (20)) up to an obvious change in notation.

Remark 5.6

In practical applications, an open challenge is to fit the model (5.14) to MRI data as in Figure 2. The statistical problem is under-specified, since there are an infinite number of choices for α(y) and m(dy) = M(y)dy that will agree with any finite data set. One reasonable approach is to fit the simplest functions α(y) and M(y) using spherical harmonics in d = 3 dimensions. For example, the data in Figure 2 can be fit using six spherical harmonics. The resulting functions α(y) and M(y) will agree exactly with the measured values of α(y) and the corresponding weights M(y) obtained from the regression lines in Figure 2, and smoothly interpolate in between.

6. Time-fractional models for DTI

Here we explore the challenge of developing effective time-fractional models for DTI. These models can be useful if the data exhibit a power law decay in S as a function of b.

Anomalous subdiffusion can be modeled using a fractional derivative in time. Given a function f(t) with Laplace transform f ~ ( s ) = ∫ 0 ∞ e - st f ( t ) dt , recall that s˜f(s) − f(0) is the Laplace transform of the first derivative f′(s). The Caputo fractional derivative ∂_t^βf(t) is defined for 0 < β < 1 as the function with Laplace transform s^β(s) − s^β−1f(0), extending the traditional form. Take Fourier and then Laplace transforms in the space-time fractional diffusion equation

∂ t β p ( x , t ) = D Δ α / 2 p ( x , t )

(6.1)

p ¯ ( k , s ) = s β - 1 s β + D ‖ k ‖ α

E β ( x ) : = ∑ n = 0 ∞ x n Γ ( 1 + β n )

E β ( x ) ~ x - 1 Γ ( 1 - β ) as x → ∞

(6.2)

Figure 3 shows a log–log plot of S versus b for the data from Figure 2, for all six directions. The straight line behavior in Figure 3 shows that a time-fractional model is a reasonable alternative to the stretched exponential, since a power law S ≈ Cb^{− β} also gives a good fit to the data. This is indicated by a straight line on the log–log plot with slope −β, since ln S ≈ ln C − β ln b. The β estimates and standard errors are listed in Table 2. OPEN IN VIEWER

OPEN IN VIEWER

Table 2.

Best-fit β values via linear regression on the data in Figure 3 for six different directions, demonstrating anisotropy.

Direction	1	2	3	4	5	6
β	0.494	0.571	0.599	0.676	0.605	0.661
Standard error	0.006	0.017	0.020	0.014	0.018	0.014

Again, it is clear that the data exhibit significant anisotropy, since the β values vary significantly with direction. For example, the value for direction 1 (corresponding to direction vector g ₁) is β = 0.494 ± 0.006 which is significantly different from the β = 0.571 value in direction 2. Since the sample size is n = 9, the 95% confidence interval is β ± 2.365 SE using the 97.5th percentile of the standard t distribution with n − 2 = 7 degrees of freedom. For example, we are 95% confident that the correct β value for direction 1 lies in the interval (0.480, 0.508).

6.1. Time-fractional Hanyga diffusion

The paper of Hanyga and Magin (2014) also proposed a time-fractional version of their anisotropic fractional diffusion model. Define the pseudo-differential operator Q on the space C₀(ℝ ^d ) of smooth functions with compact support (i.e. f(x) = 0 for all ‖x‖ ≥ M for some M > 0) such that Q f(x) has Fourier transform Q ∧ (k) f ∧ (k), where Q ∧ ( k ) is defined by (5.2). This operator can then be extended to larger spaces of functions, or even distributions (i.e. generalized functions). Since p ∧ ( k , t ) = e t Q ∧ ( k ) , it is obvious that this Fourier transform solves the ordinary differential equation

d d t p ∧ ( k , t ) = Q ∧ ( k ) p ∧ ( k , t ) ; p ∧ ( k , 0 ) ≡ 1

(6.3)

Inverting the Fourier transform shows that the probability densities p(x, t) of the stochastic process X(t) from Proposition 5.1 solve the pseudo-differential equation

∂ ∂ t p ( x , t ) = Q p ( x , t ) ; p ( x , 0 ) = δ ( x )

(6.4)

The equation (6.4) is also called a Cauchy problem (Arendt et al., 2001). In fact, if we define the semigroup

T t f ( x ) = ∫ f ( x - y ) p ( y , t ) d y

on the space L¹(ℝ ^d ) of integrable functions f: ℝ ^d  → ℝ, then Q ∧ is the generator of that semigroup (Baeumer and Meerschaert, 2001, Theorem 2.2), and T_t f(x) solves the Cauchy problem

∂ t p ( x , t ) = Q p ( x , t ) ; p ( x , 0 ) = f ( x )

(6.5)

for any f ∈ Dom( Q ∧ ), the domain of the generator. It follows from Baeumer and Meerschaert (2001, Proposition 2.1) that this semigroup T_t is strongly continuous and uniformly bounded, and that we can write the generator explicitly in the form

Q f ( x ) = ∫ y ≠ 0 ( f ( x - y ) - f ( x ) + y · ∇ f ( x ) 1 + ‖ y ‖ 2 ) φ ( d y )

for any f ∈ W²,¹(ℝ ^d ), the Sobolev space of functions in L¹(ℝ ^d ) whose first and second partial derivatives are all in L¹(ℝ ^d ).

Given any 0 < β < 1, define the Riemann–Liouville fractional derivative

D t β g ( t ) = 1 Γ ( 1 - β ) d dt ∫ 0 ∞ g ( t - s ) s - β d s

Then it follows from Baeumer and Meerschaert (2001, Theorem 3.1) that the function

q ( x , t ) = ∫ 0 ∞ p ( x , ( t / u ) β ) g β ( u ) d u

(6.6)

solves the fractional Cauchy problem

D β t q ( x , t ) = Q q ( x , t ) + t - β Γ ( 1 - β ) f ( x )

(6.7)

whenever p(x, t) solves the Cauchy problem (6.5). Here g_β(u) is the probability density function of the standard β-stable subordinator, most simply characterized in terms of its Laplace transform

∫ 0 ∞ e - st g β ( t ) dt = e - s β

(6.8)

for all s > 0, for any 0 < β < 1. A simple change of variable in the formula (6.6) reveals that

q ( x , t ) = ∫ 0 ∞ p ( x , u ) h ( u , t ) d u

(6.9)

where

h ( u , t ) = t β u - 1 - 1 / β g β ( tu - 1 / β )

(6.10)

and this leads to a stochastic solution: let D(t) be the standard β-stable subordinator, a strictly increasing infinitely divisible Lévy process such that D = D(1) has the probability density function g_β(t). Define the inverse stable process (first passage time)

E t = inf { u > 0 : D ( u ) > t }

(6.11)

and apply Corollary 3.1 from Meerschaert and Scheffler (2004) to see that the function h(u,t) in (6.10) is the probability density function of the stochastic process E_t for each t > 0. Then it follows by a standard conditioning argument that the solution q(x, t) to the fractional Cauchy problem (6.7) with the point source initial condition f(x) = δ(x) is also the probability density function of the time-changed process X(E_t), where E_t is independent from X(t). For a general initial condition f(x) that is a probability density function, the solution q(x, t) to the fractional Cauchy problem (6.7) with initial condition f(x) is the probability density function of X₀ + X(E_t), where the initial particle location X₀ has probability density function f(x). See Meerschaert and Scheffler (2008, Theorem 4.1 and Remark 4.6) for more details and extensions. Freely available R code to compute the function h(u,t) is available (Meerschaert and Sikorskii, 2012, Example 5.13) so that the solution (6.9) to the fractional Cauchy problem can be explicitly computed by numerically integrating the formula (6.9), once the probability density function p(x, t) has been computed.

The Caputo and Riemann–Liouville fractional derivatives are related by

∂ β t g ( t ) = D t β g ( t ) - t - β Γ ( 1 - β ) g ( 0 )

Then clearly one can also write the fractional Cauchy problem in a more compact form:

∂ β t q ( x , t ) = Q q ( x , t )

(6.12)

This extends the results of Hanyga (2002) for the case where β(y) is a constant.

6.2. Time-fractional Bloch–Torrey equation

Let Q ∧ be the generator of some C₀ semigroup (Arendt et al., 2001). The time-fractional Bloch–Torrey equation has been written in the literature as

∂ t β S = - ( i γ x · G ) S + D 0 Q S

(6.13)

but this form is not dimensionally correct, since the time units of ∂_t^βS are different to the units of the term ∂_tS and, more importantly, the reaction term −(iγ x  ·  G )S.

Using an idea from Baeumer et al. (2005), we can write a dimensionally correct version of the time-fractional Bloch–Torrey equation as

∂ t S = - ( i γ x · G ) S + ∇ · D ∇ ∂ t 1 - β S

Equivalently, we can write

∂ t β S = - i γ x · I t 1 - β [ G S ] + ∇ · D ∇ S

where I 1 - β is the Riemann–Liouville fractional integral defined by

I γ g ( t ) : = 1 Γ ( γ ) ∫ 0 ∞ g ( t - u ) u γ - 1 1 mu du

An alternative form is proposed by Haynga and Seredyńska (2012, equation (17)). An open challenge in the theory of DTI is to derive an analytical solution for a physically correct time-fractional Bloch–Torrey equation, suitable for clinical applications.

7. Space-variable fractional DTI models

In clinical practice, the parameters of the (fractional) Bloch–Torrey equation vary with location. Indeed, three-dimensional maps of the parameters are an important outcome of fractional DTI modeling; see for example GadElkarim et al. (2013). An open challenge is to develop the mathematical foundations of space-variable fractional DTI models. One promising approach is to use the theory of pseudo-differential operators (Schilling, 1998; Jacob, 2001). We can consider the Cauchy problem

∂ ∂ t p ( x , t ) = Q p ( x , t ) ; p ( x , 0 ) = p 0 ( x )

(7.1)

Q f ( x ) = ∫ y ≠ 0 ( f ( x - y ) - f ( x ) + y · ∇ f ( x ) 1 + ‖ y ‖ 2 ) φ ( x , d y )

Q ∧ ( x , k ) = - ∫ y ∈ Ω | y · k | α ( x , y ) m ( x , d y )

(7.2)