## Abstract

Fractional-order dynamics in physics, particularly when applied to diffusion, leads to an extension of the concept of Brownian motion through a generalization of the Gaussian probability function to what is termed anomalous diffusion. As magnetic resonance imaging is applied with increasing temporal and spatial resolution, the spin dynamics is being examined more closely; such examinations extend our knowledge of biological materials through a detailed analysis of relaxation time distribution and water diffusion heterogeneity. Here, the dynamic models become more complex as they attempt to correlate new data with a multiplicity of tissue compartments, where processes are often anisotropic. Anomalous diffusion in the human brain using fractional-order calculus has been investigated. Recently, a new diffusion model was proposed by solving the Bloch–Torrey equation using fractional-order calculus with respect to time and space. However, effective numerical methods and supporting error analyses for the fractional Bloch–Torrey equation are still limited. In this paper, the space and time fractional Bloch–Torrey equation (ST-FBTE) in both fractional Laplacian and Riesz derivative form is considered. The time and space derivatives in the ST-FBTE are replaced by the Caputo and the sequential Riesz fractional derivatives, respectively. Firstly, we derive an analytical solution for the ST-FBTE in fractional Laplacian form with initial and boundary conditions on a finite domain. Secondly, we propose an implicit numerical method (INM) for the ST-FBTE based on the Riesz form, and the stability and convergence of the INM are investigated. We prove that the INM for the ST-FBTE is unconditionally stable and convergent. Finally, we present some numerical results that support our theoretical analysis.

## 1. Introduction

Fractional calculus has been widely used in recent years in various applications in science and engineering. Fractional-order dynamics in physics, particularly when applied to diffusion, leads to an extension of the concept of Brownian motion through a generalization of the Gaussian probability function to what is termed anomalous diffusion [1–4]. A particular and very interesting class of complex phenomena arises in nuclear magnetic resonance and magnetic resonance imaging (MRI). Here, fractional calculus may help to express a physical meaning through the fractional order of the derivative operator.

As MRI is applied with increasing temporal and spatial resolution, the spin dynamics are being examined more closely; such examinations extend our knowledge of biological materials through a detailed analysis of the relaxation time distribution and water diffusion heterogeneity [1]. In many biological tissues, the diffusion-induced MR signal loss deviates from monoexponential decay, (where *D* is the diffusion coefficient and *b* is the *b*-factor), particularly at high *b* values (for example, greater than 1500 s mm^{−2} for human brain tissues) [5]. This phenomenon is sometimes referred to as anomalous diffusion, and several analytical models of non-monoexponential decay have been suggested [6,7]. A common assumption is that diffusion in two or more compartments is being observed with the compartments representing *fast* (extracellular) and *slow* (intracellular) components to the signal. The *fast* component dominates at low *b* values, and the *slow* at higher *b* values [8], leading to a biexponential form for the signal decay *S*(*b*), with a different value for the diffusion constant in each compartment
1.1where *D*_{fast}>*D*_{slow} are the diffusion constants in each compartment and *f* the volume fraction for the fast compartment.

Although intuitively appealing, there are several difficulties associated with the biexponential model. First, fitting the biexponential curve is nontrivial. It is a nonlinear fitting problem where different parameter combinations can lead to similar fits. Additionally, observed compartment sizes do not correspond to known volume fractions of intra- to extracellular space in analysed brain tissue microstructure [8,9]. Furthermore, *in vitro* experiments involving images of single oocytes also reveal non-monoexponential behaviour from the intracellular contribution alone [10], suggesting fractions obtained from biexponential fitting do not correspond to intra- and extracellular volume fractions.

Pfeuffer *et al.* [11] generalized the biexponential decay to a multi-compartmental model
1.2where *f*_{i} is the volume fraction of the *i*th compartment and the sum of *f*_{i} satisfies . The increased number of compartments provides additional degrees of freedom to fit the experimental data. However, the quality of the fit is not necessarily improved partly owing to the increased complexity in nonlinear fitting. Yablonskiy *et al.* [12] further generalized this model by replacing the discrete diffusion coefficients with a continuous distribution described by a probability function *p*(*D*). Using this model, the average cell size has been quantitatively related to measurable diffusion parameters. The exact expression of *p*(*D*), however, is unknown. The assumption of a Gaussian distribution cannot be generalized to all tissue types, considering the heterogeneous nature of tissue complexity. Instead of explicitly deriving an expression for *p*(*D*), Bennett *et al.* [13] used a stretched exponential model (equation (1.3)) to describe the diffusion-induced signal loss
1.3where DDC (coined as the distributed diffusion coefficient) is a single number representation of the diffusion coefficient distribution function *p*(*D*), and *α* is an empirical constant (0<*α*≤1). It has been demonstrated that the stretched exponential function not only fits the diffusion data from human brain tissues more precisely but also can be used to infer microscopic tissue structures through *α*, the so-called heterogeneity index [13]. The empirical stretched exponential function of Bennett *et al.* [13] was recently derived independently by Özarslan *et al.* [14] and Hall & Barrick [15], using concepts established for anomalous diffusion and fractal models.

Anomalous diffusion refers to models of diffusion in which the environment is not locally homogeneous, involving disorder that is not well approximated by assuming a uniform change in diffusion constant. Such systems include diffusion in complicated structures such as porous or fractal media. In the study performed by Hall & Barrick [16], the stretched exponential formalism was derived by recognizing that (i) the mean-squared displacement 〈*r*^{2}(*t*)〉 of diffusing molecules is related to diffusion time *t* by equation (1.4) and (ii) the dependence of apparent diffusion coefficient (ADC) on *b* can be expressed analogously to the dependence of diffusion coefficient on *t* (equation (1.5)):
1.4and
1.5where 〈*R*^{2}(*b*)〉 is the apparent mean-squared displacement, analogous to 〈*r*^{2}(*t*)〉. Equations (1.4) and (1.5) directly lead to the stretched exponential expression described by equation (1.3).

At present, a growing number of works in science and engineering deal with dynamical systems described by fractional-order equations that involve derivatives and integrals of non-integer order [16–20]. These new models are more adequate than the previously used integer-order models, because fractional-order derivatives and integrals enable the description of the memory and hereditary properties of different substances [21]. This is the most significant advantage of the fractional-order models in comparison with integer-order models, in which such effects are neglected. If the complex heterogeneity structure, such as the spatial connectivity, can facilitate movement of particles within a certain scale, fast motions may no longer obey the classical Fick's law and may indeed have a probability density function that follows a power law. Densities of *β*-stable type have been used to describe the probability distribution of these motions. The resulting governing equation of these motions is similar to the traditional diffusion equation except that the order *β* of the highest derivative is fractional. For a large number of independent solute particles, the probability propagator is replaced by the expected concentration [22]. For example, if *C*(*x*,*t*) represents the concentration of the diffusing species in one dimension, then a space and time fractional diffusion equation of the form
1.6where *K*_{x} is the generalized diffusion coefficient, is the Caputo time fractional derivative of order *α* (0<*α*<1) with respect to *t* and with the starting point at *t*=0 is defined as [21]
1.7and is the Riesz fractional derivative of order *β* (1<*β*≤2) with respect to *x*, which is defined by equation (2.1).

Although the success of using these fractional models to describe complex diffusion in biological tissues remains to be seen, they suggest a possible fractional-order dynamics in the observed diffusion-induced magnetization changes, as dictated by the Bloch–Torrey equation. Inspired by this possibility, some authors recently examined the connection between fractional-order dynamics and diffusion by solving the Bloch–Torrey equation using fractional-order calculus [1–4]. They have demonstrated that a fractional calculus-based diffusion model can be successfully applied to analysing diffusion images of human brain tissues and provide new insights into further investigations of tissue structures and the microenvironment. The following new diffusion model was proposed for solving the Bloch–Torrey equation using fractional-order calculus with respect to time and space using the Riesz formulation [1]:
1.8where *λ*=−i*γ*(**r**⋅**G**(*t*)), **G**(*t*) is the magnetic field gradient, *γ* and *D* are the gyromagnetic ratio and the diffusion coefficient, respectively, **r**=(*x*,*y*,*z*). Here, is a sequential Riesz fractional-order Laplacian operator in space [23]. *M*_{xy}(**r**,*t*)=*M*_{x}(**r**,*t*)+i*M*_{y}(**r**,*t*), where , comprises the transverse components of the magnetization; *τ*^{α−1} and *μ*^{2(β−1)} are the fractional-order time and space constants needed to preserve units (0<*α*≤1 and 1<*β*≤2). Magin *et al.* [1] have derived the analytical solutions with fractional-order dynamics in space (i.e. *α*=1, *β* an arbitrary real number, 1<*β*≤2) and time (i.e. 0<*α*<1 and *β*=2), respectively. However, effective numerical methods and supporting error analyses for the space and time fractional Bloch–Torrey equation (ST-FBTE) are still limited. This motivates us to derive an analytical solution of the fractional Laplacian formulation and an effective implicit numerical method (INM) for the ST-FBTE in Riesz form, and to study the stability and convergence of the proposed numerical method. We note in passing that the Riesz and fractional Laplacian formulations are not equivalent, except for certain special cases. For example, in one dimension with zero Dirichlet boundary conditions, these two formulations are the same [24].

In this paper, we consider the ST-FBTE with initial and boundary conditions on a finite domain.

The structure of the remainder of this paper as follows. In §2, some mathematical preliminaries are introduced. In §3, an approximate analytical solution for the ST-FBTE in fractional Laplacian form is derived. In §4, we propose an INM for the ST-FBTE in Riesz form. The stability and convergence of the INM are investigated in §§5 and 6, respectively. In §7, we exhibit the convergence rate for NME for a preliminary study based on a two-dimensional example to confirm the theoretical results reported in §6.

## 2. Preliminary knowledge

In this section, we outline important definitions and lemmas used throughout the remaining sections of this paper.

### Definition 2.1

The Riesz fractional operator **R**^{β} for *n*−1<*β*≤*n* on a infinite interval is defined as [25]
2.1where ,
2.2and
2.3

Similarly, we can define the Riesz fractional derivatives and of order *β* (1<*β*≤2) with respect to *y* and *z*.

### Definition 2.2 (Luchko & Gorenflo [26])

A real or complex-valued function *f*(*x*),*x*>0, is said to be in the space , if there exists a real number *p*>*α* such that
2.4for a function *f*_{1}(*x*) in .

### Definition 2.3 (Luchko & Gorenflo [26])

A function *f*(*x*),*x*>0, is said to be in the space , if and only if *f*^{(m)}∈*C*_{α}.

### Definition 2.4

The Mittag–Leffler function is defined as [21] and

### Lemma 2.5

*The initial value problem* (*α*>0)
2.5*where the function g*(*t*) *is assumed to lie in C*_{−1} *if* *in* *if* *and the unknown function y*(*t*) *is to be determined in the space* *has a unique solution*
2.6*where y*_{g}(*t*) *is a solution of the inhomogeneous fractional differential equation* (2.5) *with zero initial conditions and is represented in the form*
2.7*y*_{h}(*t*) *is a solution of the homogeneous fractional differential equation* (2.5) (*g*(*t*) *replaced by* 0) *with the given initial conditions, and we have*
2.8*where* .

### Proof.

See Luchko & Gorenflo [26]. □

We present our solution techniques for the ST-FBTE in the following two sections. The ST-FBTE (1.8) can be rewritten in the following form:
2.9For the numerical solution of the ST-FBTE in Riesz form, we equate real and imaginary components to express equation (2.9) as a coupled system of partial differential equations for the components *M*_{x} and *M*_{y}, namely
2.10and
2.11where *λ*_{G}=*γ*(**r**⋅**G**(*t*)).

For convenience, the ST-FBTEs (2.10) and (2.11) are decoupled, which is equivalent to solving [24] 2.12

When considering the analytical solution, the fractional Laplacian model becomes
2.13where *M* can be either *M*_{x} or *M*_{y}, and *f*(**r**,*t*)=*λ*_{G}*M*_{y}(**r**,*t*) if *M*=*M*_{x}, and *f*(**r**,*t*)=−*λ*_{G}*M*_{x}(**r**,*t*) if *M*=*M*_{y}.

## 3. Analytical solution of the ST-FBTE

We consider an analytical solution for the following ST-FBTE with initial and zero Dirichlet boundary conditions on a finite domain:
3.1
3.2
and
3.3where 0<*α*≤1, 1<*β*≤2, 0<*t*≤*T*, **r**=(*x*,*y*,*z*)∈*Ω*, *Ω* is the finite rectangular region [0,*L*_{1}]×[0,*L*_{2}]×[0,*L*_{3}] and *Γ* is the boundary of *Ω*.

Following the work of Ilic *et al.* [27] and Yang *et al.* [24], we set
3.4and
3.5where , and are the eigenvalues and corresponding eigenfunctions of the three-dimensional Laplacian (−*Δ*) for *n*,*m*,*l*=1,2,… and

Substituting (3.4) and (3.5) into (3.1), we have 3.6namely, 3.7

Since *M*(**r**,*t*) must also satisfy the initial condition (3.2),
3.8we obtain
3.9

For each value of *n*, *m* and *l*, (3.7)–(3.9) comprise a three-dimensional fractional initial value problem. According to lemma 2.5 and definition 2.4, the three-dimensional fractional initial value problem has the solution
3.10in which
3.11and
3.12

Hence, the analytical solution of equations (3.1)–(3.3) is given by 3.13

## 4. Implicit numerical method for the ST-FBTE

We now propose an INM for the following ST-FBTE in Riesz form with initial and boundary conditions on a finite domain:
4.1
4.2
and
4.3where 0<*α*≤1, 1<*β*≤2, 0<*t*≤*T*, **r**=(*x*,*y*,*z*)∈*Ω*, *Ω* is the finite rectangular region [0,*L*_{1}]×[0,*L*_{2}]×[0,*L*_{3}] and *Γ* is the boundary of *Ω*.

Let *h*_{x}=*L*_{1}/*N*_{1}, *h*_{y}=*L*_{2}/*N*_{2}, *h*_{z}=*L*_{3}/*N*_{3} and *τ*=*T*/*N* be the spatial and time steps, respectively. At a point (*x*_{i},*y*_{j},*z*_{k}) at the moment of time *t*_{n} for and , we denote the exact and numerical solutions *M*(**r**,*t*) as *u*(*x*_{i},*y*_{j},*z*_{k},*t*_{n}) and , respectively.

Firstly, adopting the discrete scheme in Shen *et al.* [28], we discretize the Caputo time fractional derivative of *u*(*x*_{i},*y*_{j},*z*_{k},*t*_{n+1}) as
4.4where *b*_{l}=(*l*+1)^{1−α}−*l*^{1−α}, *l*=0,1,…,*N*.

Using the relationship between the Riemann–Liouville derivative and the Grünwald–Letnikov scheme, we discretize the Riesz fractional derivative by the shifted Grünwald–Letnikov scheme in Yang *et al.* [24]
4.5and
4.6where the coefficients are defined by
4.7Similarly,
4.8
4.9
4.10
and
4.11

Thus, we can derive the implicit numerical scheme: 4.12

We then have the following implicit difference approximation:
4.13with
4.14and
4.15where *μ*_{0}=*τ*^{α}*Γ*(2−*α*)/*K*_{α}, , , , and noting that coefficients *μ*_{0}>0, *μ*_{1},*μ*_{2},*μ*_{3}<0 for 0<*α*≤1 and 1<*β*≤2.

### Lemma 4.1 (Liu *et al*. [29])

*The coefficients b*_{l}, *l*=0,1,2,… *satisfy*

(i)

*b*_{0}=1,*b*_{l}>0*for l*=1,2,…(ii)

*b*_{l}>*b*_{l+1}*for l*=0,1,2,…

## 5. Stability of the implicit numerical method for the ST-FBTE

We prove the stability of the INM for the ST-FBTE.

Let be the approximate solution of the INM (4.13)–(4.15), , and .

Assuming and using mathematical induction, we have the following theorem.

### Theorem 5.1

*The INM defined by (4.13)–(4.15) is unconditionally stable, and*
5.1

### Proof.

According to (4.13)–(4.15), the error satisfies 5.2

When *n*=0, assume that . Using lemmas 4.1 and 4.2, and noting that *μ*_{1},*μ*_{2},*μ*_{3}<0, we have

Now, suppose that By assuming that , using lemmas 4.1 and 4.2 again, we obtain

Hence the INM defined by (4.13)–(4.15) is unconditionally stable. □

## 6. Convergence of the implicit numerical method for the ST-FBTE

We prove the convergence of the INM for the ST-FBTE.

Setting , and we denote , then **R**^{0}=**0**. Here, **R**^{n} and **0** are ((*N*_{1}−1)×(*N*_{2}−1)×(*N*_{3}−1)) vectors, respectively.

From (4.1)–(4.15), the error satisfies
6.1for *i*=1,2,…,*N*_{1}−1, *j*=1,2,…,*N*_{2}−1, *k*=1,2,…,*N*_{3}−1.

Assuming , and using mathematical induction, we have the following theorem.

### Theorem 6.1

*The implicit difference approximation defined by (4.13)–(4.15) is convergent, and there is a positive constant C***, such that
*6.2

### Proof.

When *n*=0, assume that . Using equation (6.1), lemmas 4.1 and 4.2, we have

Now, suppose that , *m*=1,2,…,*n*. By assuming that . Using lemmas 4.1 and 4.2 and equation (6.1) again, we have
We note that , there exists a positive constant *C*_{2}, such that .

Finally, note that *nτ*≤*T* is finite, so there exists a positive constant *C**, such that for *n*=0,1,2,….

## 7. Numerical results

Owing to the complexity of the analytical solution given in this paper, and the computational overheads necessary to perform the simulations for NME in three dimensions, we present here a preliminary study based on a two-dimensional example to confirm our theoretical analysis.

### Example 7.1

The following ST-FBTE with initial and boundary conditions on a finite domain is considered:
7.1
7.2
and
7.3where
7.4and 0<*α*≤1, 1<*β*≤2, *t*>0, **r**=(*x*,*y*)∈*Ω*, *Ω* is the finite rectangular region [0,1]×[0,1] and *Γ* is the boundary of *Ω*.

The exact solution of this problem is *M*(**r**,*t*)=*t*^{α+β}*x*^{2}(1−*x*)^{2}*y*^{2}(1−*y*)^{2}, which can be verified by substituting directly into (7.1).

The maximum absolute error between the exact solution and the numerical solutions by INM, with spatial and temporal steps at time *t*=1.0 when *K*_{α}=1.0,*K*_{β}=0.5,*α*=0.8,*β*=1.8, is listed in table 1. Table 2 shows the maximum absolute error between the exact solution and the numerical solutions by INM, with spatial steps at time *t*=1.0 when *K*_{α}=1.0,*K*_{β}=0.5,*α*=0.8,*β*=1.8.

From table 1, it can be seen that From table 2, it can be seen that

This is in good agreement with our theoretical analysis, namely that the convergence order of the numerical method INM for this problem is *O*(*τ*^{2−α}+*h*_{x}+*h*_{y}).

We now exhibit in example 7.2 the result of INM for the ST-FBTE with a nonlinear source term.

### Example 7.2

ST-FBTE with initial and boundary conditions on a finite domain:
7.5
7.6
and
7.7where the nonlinear source term is Fisher's growth equation *f*(*M*,**r**,*t*)=0.25*M*(**r**,*t*)[1−*M*(**r**,*t*)], and 0<*α*≤1, 1<*β*≤2, *t*>0, **r**=(*x*,*y*)∈*Ω*, *Ω* is the finite rectangular region [0,1]×[0,1] and *Γ* is the boundary of *Ω*.

The solution profiles of (7.5) by INM, with spatial and temporal steps , at time with *K*_{α}=1.0,*t*_{final}=1.0 for different *α*,*β* and *K*_{β} are shown in figure 1. From figure 1, it can be seen that the coefficient *K*_{β} impacts on the solution profiles of (7.5), whereby a larger value of *K*_{β} produces more diffuse profiles. In figure 2, we illustrate the effect of the fractional order in space for this problem, with spatial and temporal steps , at time with *K*_{α}=1.0,*K*_{β}=1.0,*t*_{final}=1.0 for *β* fixed at 2 and *α* varying. Here, we can see that reducing the value of *α* leads to a much sharper central peak. In figure 3, we illustrate the effect of the fractional order in time for this problem, with spatial and temporal steps , at time with *K*_{α}=1.0,*K*_{β}=1.0,*t*_{final}=1.0 for *α* fixed at 1 and *β* varying. Again we see the peak becoming sharper as *β* is decreased. We now see that a much sharper profile is obtained when we vary *α* with *β* fixed than when we vary *β* with *α* fixed. Finally, the effect of the fractional order in both time and space is illustrated in figure 4.

Overall, the coefficient *K*_{β} has a significant effect, especially for fractional values of *α*. The most significant effects in terms of the spikiness of the profile occur as both *α* and *β* are simultaneously reduced.

## 8. Conclusions

In this paper, an analytical solution and an effective INM for solving the fractional Bloch–Torrey equation in fractional Laplacian and Riesz forms, respectively, have been derived. The stability and convergence of the INM are analysed systematically. We have used our numerical method to simulate a problem of practical importance involving a nonlinear source term. Our results have highlighted the impact of the fractional indices on the shape of the solution profile. Finally, we note, however, that the complexity of our solutions and the computational overheads of our methods are such that we plan in the future to investigate alternative solution strategies for solving ST-FBTEs such as those based on alternating direction implicit methods.

## Acknowledgements

The authors gratefully acknowledge the help in their work by Prof. Kerrie Mengersen from Queensland University of Technology. Q.Y. also acknowledges the Centre for Complex Dynamic Systems Control for offering financial support for his PhD scholarship.

## Footnotes

One contribution of 14 to a Theme Issue ‘Fractional calculus and its applications’.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.