- Email: [email protected]
Fractional Bloch equation with delay
Computers and Mathematics with Applications 61 (2011) 1355–1365
Contents lists available at ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
Fractional Bloch equation with delay Sachin Bhalekar a , Varsha Daftardar-Gejji a , Dumitru Baleanu b,c , Richard Magin d,∗ a
Department of Mathematics, University of Pune, Pune - 411007, India
b
Department of Mathematics and Computer Science, Faculty of Arts and Sciences, Cankaya University - 06530, Ankara, Turkey
c
Institute of Space Sciences, P.O. Box, MG-23, R 76900, Magurele-Bucharest, Romania
d
Department of Bioengineering, University of Illinois, 851 S. Morgan St., Chicago, 60607, USA
article
abstract
info
Article history: Received 10 July 2010 Received in revised form 30 December 2010 Accepted 31 December 2010
In this paper we investigate a fractional generalization of the Bloch equation that includes both fractional derivatives and time delays. The appearance of the fractional derivative on the left side of the Bloch equation encodes a degree of system memory in the dynamic model for magnetization. The introduction of a time delay on the right side of the equation balances the equation by also adding a degree of system memory on the right side of the equation. The analysis of this system shows different stability behavior for the T1 and the T2 relaxation processes. The T1 decay is stable for the range of delays tested (1–100 µs), while the T2 relaxation in this model exhibited a critical delay (typically 6 µs) above which the system was unstable. Delays are expected to appear in NMR systems, in both the system model and in the signal excitation and detection processes. Therefore, by including both the fractional derivative and finite time delays in the Bloch equation, we believe that we have established a more complete and more realistic model for NMR resonance and relaxation. © 2011 Elsevier Ltd. All rights reserved.
Keywords: Fractional calculus Bloch equation Delay
1. Introduction The classical description of nuclear magnetic resonance (NMR) – the phenomena underlying magnetic resonance imaging (MRI) – is summarized in vector form by the Bloch equation [1,2], →
dM dt
→
→
=γ M× B −
(Mz − M0 ) T1
iz −
(Mxix + Myiy ) T2
,
→
(1.1) →
where M (Mx , My , Mz ) represents the time-varying magnetization (Amps/s), M0 the equilibrium magnetization, B (Bx , By , Bz ) the applied radio frequency (RF), gradient and static magnetic fields (Tesla), γ is the gyromagnetic ratio (42.57 MHz/T, for spin 1/2 protons), and T1 , T2 are the spin–lattice and the spin–spin relaxation times, respectively. For homogeneous (e.g., over a 1–2 mm3 voxel), and isotropic materials with a single spin component (typically water protons), the Bloch equation prescribes the dynamic balance between externally applied magnetic fields and internal sample relaxation times. This dynamic balance is the basis for pulse sequence design, signal acquisition, image reconstruction and, in the case of MRI, tissue contrast [3–5]. Fractional order generalization of the Bloch equation has been undertaken by several groups [6–13] to account for the anomalous relaxation and anomalous diffusion observed in NMR studies of complex materials — typically gels, emulsions,
∗
Corresponding author. Tel.: +1 312 413 5528; fax: +1 312 996 5921. E-mail addresses: [email protected] (S. Bhalekar), [email protected], [email protected] (V. Daftardar-Gejji), [email protected] (D. Baleanu), [email protected] (R. Magin). 0898-1221/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.camwa.2010.12.079
1356
S. Bhalekar et al. / Computers and Mathematics with Applications 61 (2011) 1355–1365
porous composites and biological tissues. A common feature of such complex materials is a ‘‘mesoscopic’’ structure intermediate in scale between the molecular and the macroscopic regimes. NMR spectroscopy and MRI are powerful experimental tools for probing the organization (ordered/disordered) and the dynamics (phase transitions, diffusion, permeability) of mesoscopic structures. High resolution MR microscopy, for example, provides high contrast images of phase separation, particle aggregation, and macromolecular coalescence, while NMR diffusion studies give a direct measure of material porosity and tortuosity. However, the application of these experimental tools and the analysis of the acquired data are highly dependent upon the theoretical assumptions underlying the Bloch equation: typically Gaussian spin dynamics, Brownian particle motion, and first order exponential relaxation. To the extent that these assumptions apply to a particular material or to a transition between different material phases, conventional NMR analysis is valid and appropriate. Nevertheless, accumulating experimental evidence on complex materials, acquired using dielectric spectroscopy, dynamic light scattering as well as ultrasonic and viscoelastic measurements, strongly suggests anomalous dynamic behavior. This anomalous behavior appears to reflect distributions of relaxation times, and multi-scale phenomena that in some cases suggest a fractallike structure, nonlocal interactions, fading memory and aging. Hence, it is anticipated that NMR and MRI measurements of complex materials with mesoscopic states or configurations will show similar – fractional order – dynamic behavior. The largely unresolved question concerning the fractional or ‘‘generalized’’ Bloch equation is how to modify it to relate the fractional order of the operators to the spin dynamics at the molecular, mesoscopic, and macroscopic scales. NMR is a widely used experimental technique because it effectively probes time scales from the nanosecond (spin precession) to the microsecond (molecular reorientation) to the millisecond (spin–spin dephasing) and up to the second (spin–lattice relaxation) and beyond (diffusion). Using images or localized spectroscopy, NMR analysis assigns to each region of interest (pixel or voxel) specific experimental values (e.g., diffusion coefficient, spin–spin and spin–lattice relaxation times, correlation coefficient, chemical shift). This process works well for imaging macroscopic regions (i.e., voxels on the order of 1 mm3 ) or for spectroscopy on pure samples in uniform magnetic fields (i.e., one and two dimensional NMR spectra). However, pushing the imaging resolution into the microscopic or mesoscopic regimes or pulling complex spectra from macroscopic mixtures or tissues exposes fundamental limitations in either signal-to-noise (SNR) or spectral resolution. These problems, unfortunately, cannot always be overcome by taking more signal averages or by the use of extremely high field NMR systems. The goal of fractional order NMR modelling is to extend the experimental window of MR techniques in both space and time, not by simply increasing system resolution or bandwidth, but by appending to the governing equations a fractional order interpretation for the observed phenomena. The power of fractional calculus is that it conveniently assimilates into the order of the operator a fractal network of lacunae extending over a wide distribution of relaxation times or diffusion coefficients. Such conciseness comes at the cost of relinquishing the specific sub-compartment values for individual relaxation times or diffusion constants, but provides the benefit of capturing hierarchical, multi-scale phenomena in space and time. In previous work we have investigated spatial complexity via fractional order space derivatives [11] and anomalous temporal relaxation through fractional order time derivatives [12]. Both approaches appear promising for evaluating porous, layered or semi-permeable composite materials. In this paper we investigate – for the case of relaxation – the effect of introducing a short time delay into the fractional order model. Such time delays may perhaps be used to capture some of the longer time dynamic changes in complex systems associated with aging in both the experimental and system senses. In this paper we investigate – for the case of relaxation – the effect of introducing a short time delay into the fractional order model. Such time delays are common in NMR pulse programming where they are used to generate coherence in the detected signals, and to account for imperfections in the applied gradients and radio frequency pulses. Such time delays may also capture some of the longer time dynamic changes in complex systems associated with aging in both the experimental and system senses. Fractional order operators in time are well known to incorporate system memory; this arises through the inclusion of a power law function in the time domain convolution operator definition of fractional order integration and differentiation. Fractional order operators in space – in a similar manner – extend spatial derivatives through a non-local or distributed spatial gradient that captures a wider range of sample structure in its formal definition. The delay parameter investigated in this paper modifies the time window for the fractional order time derivative by allowing not only fading memory of an earlier state or time, but by also introducing specific information about the initial conditions at a time specified in the lower limit of the convolution operator. The rationale for this approach is that we hope the generalization will extend the modelling capabilities of the Bloch equation so that transient mesoscopic structures in complex materials can be better visualized. The organization of this paper is as follows: in Section 2 (Preliminaries) we present the fundamental definitions for the fractional order operators to be used; in Section 3 (Methods) we outline the numerical method used to solve the fractional Bloch equation with delay; in Section 4 (Theory) we describe the fractional Bloch equation with both fading memory and delay included; in Section 5 (Results) we present our findings for examples selected to span a range of fractional orders and short time delay; and in the last two sections we summarize our contribution by comparing and contrasting our results with our previous work and the work of others. 2. Preliminaries Basic definitions and properties of fractional derivative/integrals are given below [14–16].
S. Bhalekar et al. / Computers and Mathematics with Applications 61 (2011) 1355–1365
1357
Definition 2.1. A real function f (t ), t > 0 is said to be in space Cα , α ∈ ℜ if there exists a real number p (>α), such that f (t ) = t p f1 (t ) where f1 (t ) ∈ C [0, ∞). Definition 2.2. A real function f (t ), t > 0 is said to be in space Cαm , m ∈ N
{0} if f (m) ∈ Cα .
Definition 2.3. Let f ∈ Cα and α ≥ −1, then the (left-sided) Riemann–Liouville integral of order µ, µ > 0 is given by µ
I f (t ) =
t
∫
1
Γ (µ)
(t − τ )µ−1 f (τ )dτ ,
t > 0.
(2.1)
0
m Definition 2.4. The (left sided) Caputo fractional derivative of f , f ∈ C− 1, m ∈ N
Dµ f (t ) =
dm dt m
f (t ),
= I m−µ
{0}, is defined as:
µ=m
dm f ( t ) dt m
,
m − 1 < µ < m, m ∈ N .
(2.2)
Note that for m − 1 < µ ≤ m, m ∈ N, I µ Dµ f (t ) = f (t ) −
m−1
− dk f k =0
I µt ν =
tk
(0) , k!
dt k
(2.3)
Γ (ν + 1) t µ+ν . Γ (µ + ν + 1)
(2.4)
Stability analysis Deng et al. [17] have studied the stability of n-dimensional linear fractional differential equation with time delays for the following system: Dα1 x1 (t ) = a11 x1 (t − τ11 ) + a12 x2 (t − τ12 ) + · · · + a1n xn (t − τ1n )
Dα2 x2 (t ) = a21 x1 (t − τ21 ) + a22 x2 (t − τ22 ) + · · · + a2n xn (t − τ2n )
··· Dαn xn (t ) = an1 x1 (t − τn1 ) + an2 x2 (t − τn2 ) + · · · + ann xn (t − τnn ),
(2.5)
where αi ∈ (0, 1), xi (t ), xi (t − τij ) state variables and τij are delay terms. It can be shown that if all the roots of the characteristic equation det(∆(λ)) = 0, where
α λ 1 − a11 e−λτ11 −a21 e−λτ21 ∆(λ) = . ..
−a12 e−λτ12 λα2 − a22 e−λτ22 .. .
−an1 e−λτn1
−an2 e−λτn2
··· ··· .. . ···
−a1n e−λτ1n −a2n e−λτ2n .. .
λαn − ann e−λτnn
,
(2.6)
have negative real parts, then the system (2.5) is asymptotically stable. 3. Algorithm for solving delay differential equations of fractional order Due to the nonlocal nature of fractional derivative operator, the numerical methods for solving ordinary differential equations have to be modified to cater for solving fractional differential equations (FDEs) [18–20]. Bhalekar and DaftardarGejji have modified Adams–Bashforth method to solve delay differential equations of fractional order (FDDE) [21]. The method is described below. Consider the FDDE Dα y(t ) = f (t , y(t − τ )), y(t ) = g (t ),
t ∈ [0, T ], 0 < α ≤ 1,
t ∈ [−τ , 0],
(3.1) (3.2)
where Dα denotes the Caputo fractional derivative of order α ∈ (0, 1]. Applying fractional integration Itα on both sides of (3.1) and using (3.2) we obtain y(t ) = g (0) +
1
Γ (α)
t
∫ 0
(t − ξ )α−1 f (ξ , y(ξ − τ ))dξ .
(3.3)
1358
S. Bhalekar et al. / Computers and Mathematics with Applications 61 (2011) 1355–1365
Consider a uniform grid {tn = nh : n = −k, −k + 1, . . . , −1, 0, 1, . . . , N } where k and N are integers such that h = T /N and h = τ /k. Let yh (tj ) = g (tj ),
j = −k, −k + 1, . . . , −1, 0
(3.4)
and note that yh (tj − τ ) = yh (jh − kh) = yh (tj−k ),
j = 0, 1, . . . , N .
(3.5)
Suppose we have already calculated approximations yh (tj ) ≈ y(tj ), (j = −k, −k + 1, . . . , −1, 0, 1, . . . , n) and we want to calculate yh (tn+1 ) using tn+1
∫
1
(tn+1 − ξ )α−1 f (ξ , y (ξ − τ ))dξ . (3.6) Γ (α) 0 We use approximations yh (tn ) for y(tn ) in (3.6). The integral on right hand side of (3.6) is evaluated using the product y(tn+1 ) = g (0) +
trapezoidal quadrature formula. We obtain yh (tn+1 ) = g (0) +
= g (0) +
hα
Γ (α + 2) hα
Γ (α + 2)
f (tn+1 , yh (tn+1 − τ )) +
f (tn+1 , yh (tn+1−k )) +
hα
n −
Γ (α + 2)
j =0
hα
n −
Γ (α + 2)
j =0
aj,n+1 f (tj , yh (tj − τ ))
aj,n+1 f (tj , yh (tj−k )),
(3.7)
where aj,n+1 are given by aj,n+1
nα+1 − (n − α)(n + 1)α , = (n − j + 2)α+1 + (n − j)α+1 − 2(n − j + 1)α+1 , 1,
if j = 0, if 1 ≤ j ≤ n, if j = n + 1.
(3.8)
Note that, unlike in the non-delayed case the right hand side of the Eq. (3.7) does not contain the unknown term yh (tn+1 ), so it is not needed to calculate the predictor term. 4. Fractional Bloch equation with delay We apply the method described in Section 3 for solving fractional order Bloch equation involving time delay T −α Dα Mx (t ) = ω˜0 My (t − τ ) −
Mx (t − τ ) T2
T −α Dα My (t ) = −ω˜0 Mx (t − τ ) − T −α Dα Mz (t ) = Mx (t ) = 0, where ω˜0 =
ω0
T α−1
M0 − Mz (t − τ ) T1 My (t ) = 100,
= (ω0 T )T −α ,
1 T′
=
1
,
My (t − τ ) T2
(4.1)
,
(4.2)
,
(4.3) Mz (t ) = 0,
1 T α−1 T
= 1
for t ≤ 0,
T −α 1 T , T′ T1 2
=
(4.4)
1 T α−1 T
2
=
T −α T . We assume that T T2 1 1 0 T ′ and T ′ have the unit
ω˜0 = 3 × π , T1′ = 1, T2′ = 20 × 10−3 , M0 = 100. We mention that ω˜ ,
1
2
= 1.0. In our case we have of (s)−α .
We use the stability analysis given in Section 2 to study the behavior of the system (4.1)–(4.3). Note that, Eq. (4.3) can be transformed to Dα Mz1 (t ) = −Mz1 (t − τ ),
(4.5)
where Mz1 (t ) = (Mz (t ) − M0 )/T1 . The characteristic equation for this system becomes:
λ2α +
2 α −λτ λ e + T2
1
1
e−2λτ
(λα + e−λτ ) = 0. (4.6) ω02 T22 + 100λα e−λτ + (102400π 2 + 2500)e−2λτ ) in (4.6) corresponds to Eq. (4.1)–(4.2) and the factor +
Note that the factor (λ2α (λα + e−λτ ) corresponds to a decoupled Eq. (4.5). We rename (4.1)–(4.2) as subsystem (A) and (4.5) as subsystem (B). The system changes its behavior (from damped oscillations to oscillations with an increasing height) for a suitable value of delay. Note that when the delay is zero, the solutions to this system of fractional order differential equations are given in [12]. When α = 1 and the delay is zero, we find the classical results. 5. Results The observations are summarized in Figs. 1–6. We consider three values of α for both system A and system B. In this example we assume a perfect spin excitation ( π2 )x so that Mx (0) = 0, My (0) = M0 , Mz (0) = 0. This system should return to equilibrium at long times compared with T1 and T2 , that is Mx (∞) = 0, My (∞) = 0, Mz (∞) = M0 .
S. Bhalekar et al. / Computers and Mathematics with Applications 61 (2011) 1355–1365
a
b
Mx
Mz
c
My 100
80
50
60
100 50 t 0.01
0.02
0.03
0.01
0.04
0.03
0.04
40
t
20
-50
-50
-100
-100
t 1
My
d
0.02
1359
2
3
4
e
100
50
-100
-50
50
100
Mx
100 50
M z 642 0 -100
-50
0 M y -50
-50
0 -100
Mx
50
-100
100
Fig. 1. A set of plots of magnetization versus time for the Mx , My and Mz relaxation dynamics following a 90x degree pulse excitation (a, b, and c respectively) and a two-dimensional plot of Mx versus My (d) and a three-dimensional plot of Mx , My and Mz (e). In this figure the time delay was 0.00006 s and the value of α was 0.8.
a
b
Mx
150000
150000
100000
100000
50000
50000 0.01
0.02
0.03
0.04
c
-100000
Mz 80
t
-50000
0.01
0.02
0.03
0.04
t
60
-50000
40
-100000
20
-150000
-150000
d
My
t 1
My
2
3
4
e
150000 100000 50000 Mx
-150000 -100000 -50000
50000 100000 150000 2000
-50000 -100000 -150000
Mz 6402 -20000
0 My 0 Mx
-20000 20000
Fig. 2. A set of plots of magnetization versus time for the Mx , My and Mz relaxation dynamics following a 90x degree pulse excitation (a, b, and c respectively) and a two-dimensional plot of Mx versus My (d) and a three-dimensional plot of Mx , My and Mz (e). In this figure the time delay was 0.00007 s and the value of α was 0.8.
In the case of α = 0.80: for τ = 0.00006, roots of the characteristic Eq. (4.6) corresponding to system (A) are −143.972 ± 5729.25i and corresponding to system (B) are −163062 ± 123853i. In view of the analysis given in Section 2,
1360
S. Bhalekar et al. / Computers and Mathematics with Applications 61 (2011) 1355–1365
a
c
b
Mx
100
100
80
50
60
50 0.01
0.02
0.03
t 0.04
0.01
-50
-50
-100
-100
My
d
Mz
My
0.02
40
t 0.04
0.03
20 t 1
2
3
4
e
100
50
-100
-50
50
Mx
100
100 50
M z 43210
-50
0 M y
-100 -50
-50
0 -100
50
Mx
-100 100
Fig. 3. A set of plots of magnetization versus time for the Mx , My and Mz relaxation dynamics following a 90x degree pulse excitation (a, b, and c respectively) and a two-dimensional plot of Mx versus My (d) and a three-dimensional plot of Mx , My and Mz (e). In this figure the time delay was 0.00009 s and the value of α was 0.9.
a
b
Mx
c
My
200
200
100
100
Mz 80
0.01
0.02
0.03
0.04
60
t
0.01
-100
-100
-200
-200
0.02
0.03
0.04
t
40 20 t 1
2
3
4
My
d
e
200
100
-200
-100
100
-100
200
Mx 200 100
Mz 43210 -200
0 My
-100 Mx -200
-100
0 100 200
-200
Fig. 4. A set of plots of magnetization versus time for the Mx , My and Mz relaxation dynamics following a 90x degree pulse excitation (a, b, and c respectively) and a two-dimensional plot of Mx versus My (d) and a three-dimensional plot of Mx , My and Mz (e). In this figure the time delay was 0.00010 s and the value of α was 0.9.
the system has stable solutions in this case. The numerical results also have the same conclusions. The damped oscillations of Mx and My are shown in Fig. 1(a) and (b), whereas Fig. 1(c) shows the stable solution for Mz .
S. Bhalekar et al. / Computers and Mathematics with Applications 61 (2011) 1355–1365
1361
c a
b
Mx
Mz 100
My
100
100
80
50
50
60
0.01
0.02
0.03
t
0.04
0.01
-50
-50
-100
-100
0.02
0.03
0.04
40
t
20 1
2
3
t
4
My
d
e
100
50
-50
50
100
Mx 100 50
Mz 3210
-50
0 My
-50 Mx
-50
0 50
100 -100
-100
Fig. 5. A set of plots of magnetization versus time for the Mx , My and Mz relaxation dynamics following a 90x degree pulse excitation (a, b, and c respectively) and a two-dimensional plot of Mx versus My (d) and a three-dimensional plot of Mx , My and Mz (e). In this figure the time delay was 0.00004 s and the value of α was 1.0.
b
a
Mx
c
My 150
150
80
100
100
60
50
50 0.01
0.02
0.03
0.04
t
0.01
-50
-50
-100
-100
-150
-150
d
Mz 100
0.02
0.03
0.04
t
40 20 1
2
3
4
t
e
My 150 100 50
-150
-100
-50
50
100
150
Mx
-50
10 0 Mz 43210
-100 -150
0 My
-100 Mx
0
-100 100
Fig. 6. A set of plots of magnetization versus time for the Mx , My and Mz relaxation dynamics following a 90x degree pulse excitation (a, b, and c respectively) and a two-dimensional plot of Mx versus My (d) and a three-dimensional plot of Mx , My and Mz (e). In this figure the time delay was 0.00006 s and the value of α was 1.0.
1362
S. Bhalekar et al. / Computers and Mathematics with Applications 61 (2011) 1355–1365
When τ is increased to 0.00007, the roots of characteristic equation become 184.389 ± 5575.72i (system (A)) and −135039 ± 9804.27i (system (B)). Since the real parts of roots for system (A) are non-negative, we cannot use the theory in Section 2. The numerical results shows that the system (A) becomes unstable for these values of parameters. In Fig. 2(a) and (b) the oscillations of Mx and My are with increasing height. This behavior is not normally observed in NMR experiments without delay unless addition RF pulses are applied that generate so-called spin echoes. The system (B) is stable (Fig. 2(c)) in this case. In the case of α = 0.90: damped oscillations of Mx and My are observed for τ = 0.00009 (Fig. 3(a), (b)). This observation by numerical methods is supported by the theory also, as the roots of the characteristic equation corresponding to system (A) are −26.4959 ± 2175.72i. Also for system (B), the roots are −116677 ± 3817.75i and the system is stable (Fig. 3(c)). For τ = 0.00010, the characteristic roots become 23.2458 ± 2164.4i (system (A)) and −103972 ± 3439.19i (system (B)). The numerical solutions (Fig. 4(a),(b)) shows that the system (A) is unstable and system (B) (Fig. 4(c)) is stable. In the case of α = 1.00: it can be observed from analytically and numerically that the system is stable for τ = 0.00004 (Fig. 5(a)–(c)). This is the expected classical behavior. For τ = 0.00006, oscillations of Mx , My with increasing height are observed (Fig. 6(a), (b)) and Mz is stable (Fig. 5(d)). Thus, even for the integer order time derivative, a short time delay is sufficient to bring about a shift in the dynamical behavior of the system. In addition all three values of α(0.80, 0.90, 1.0) a very short time delay (τ = 0.000001 s) yielded only stable anomalous relaxations via Mittag-Leffler functions. In this example we have assumed ω0 = 160 × 2π rad/s. This corresponds to a sinusoidal oscillator with a period of (T0 = 2ωπ ) of 6.25 ms, while the delay typically ranges from 0 40–100 µs. No available behavior was observed for the delay of 1 µs. 5.1. Comparison with existing method In this subsection, we compare our results with the existing discretization method discussed in [15,22,23] by Podlubny and co-workers. Let f (t ) be a function defined on [0, b]. Consider the equidistant nodes ti = ih (i = 0, 1, . . . , N ) with the step h, where t0 = 0 and tN = b. If f (j) (0) = 0, j = 0, 1, . . . , m then we can approximate fractional order derivative Dα f (t ) of order α (m − 1 < α ≤ m) at point ti as follows: Dα f (t ) |t =ti ≈ h−α
i α − (−1)r fi−r ,
(5.1)
r
r =0
where
α r
=
Γ (α + 1) r !Γ (α − r + 1)
and fi = f (ti ). Using the transformation x(t ) = Mx (t ), y(t ) = My (t ) − 100, z (t ) = Mz (t ), we can write the system (4.1)–(4.3) as Dα x(t ) = ω˜0 (y(t − τ ) + 100) −
M0 T1
x(t ) = 0,
−
z (t − τ ) T1
T2
y(t − τ ) + 100
Dα y(t ) = −ω˜0 x(t − τ ) − Dα z (t ) =
x(t − τ )
T2
,
(5.2)
,
(5.3)
,
y(t ) = 0,
(5.4) z (t ) = 0,
for t ≤ 0, 0 < α ≤ 1.
(5.5)
Consider a uniform grid {tn = nh : n = −k, −k + 1, . . . , −1, 0, 1, . . . , N } where k and N are integers such that h = b/N and h = τ /k. Using (5.1) we can write the approximation of system (5.2)–(5.5) as h−α
i −
(−1)r
r =0
h−α
i −
(−1)r
r =0
h−α
i − r =0
(−1)r
α r
α r
α r
xi−r = ω˜0 (yi−k + 100) − yi−r = −ω˜0 xi−k −
zi−r =
M0 T1
−
zi−k T1
,
x i −k T2
yi−k + 100 T2
, ,
xi = yi = zi = 0, for i ≤ 0.
S. Bhalekar et al. / Computers and Mathematics with Applications 61 (2011) 1355–1365
a
b
Mx
1363
My
100 100 50 50
0.01
0.02
0.03
0.04
0.01
t
0.02
0.03
0.04
t
-50 -50 -100 -100
c
d
Mz
My 100
80 50 60
-100
40
-50
50
100
Mx
-50 20 -100 t 1
2
3
4
Fig. 7. Plots of magnetization versus time for the Mx , My and Mz relaxation dynamics following a 90x degree pulse excitation (a, b, and c respectively) and a two-dimensional plot of Mx versus My (d), with time delay 0.00006 s and α = 0.8 using the method given in Section 5.1.
Hence we can derive the following algorithm for obtaining a numerical solution xi = yi = zi = 0,
for i ≤ 0,
[
xi−k xi = hα ω˜0 (yi−k + 100) − T2
[
yi−k + 100 yi = hα −ω˜0 xi−k − T2 z i = hα
[
M0 T1
−
zi−k T1
] −
] −
i −
(−1)r
r =1
] −
r
xi−r ,
i α − (−1)r yi−r , r =1
i α − (−1)r zi−r , r =1
α
r
r
(5.6)
where i = 1, 2, . . . , N. In Fig. 7, we have plotted the case α = 0.8, τ = 0.0006 using the algorithm (5.6). Comparing these figures with Fig. 1, we can conclude that the results by our method agree well with the existing discretization method. From this we can conclude that our results match well with the classical method and are consistent with the stability analysis discussed in Section 2. 6. Discussion In this paper we have investigated a generalization of the Bloch equation that includes both fractional derivatives and time delays. The combination of these two powerful techniques is justified for by both theoretical and experimental reasons. First, in theory, the fractional derivative appears on the left side of the Bloch equation, hence the equation is asymmetric with the memory encoded in the definition of the fractional derivative contributing only to one part of the dynamic model for magnetization. In order to account for this discrepancy, we have introduced system memory on the right side of the Bloch equation through a finite time delay — tau. Second, experimentally, delays are expected to appear in NMR systems, delay in both the system model and in the signal excitation and detection processes. Therefore, by including both the fractional derivative and finite time delays in the Bloch equation, we believe that we have established a more complete and more realistic model for NMR resonance and relaxation. In the following we will describe more completely our rationale for this generalization and the implications of our results for the analysis of NMR signals. The application of fractional
1364
S. Bhalekar et al. / Computers and Mathematics with Applications 61 (2011) 1355–1365
calculus to NMR is relatively recent and it can be seen in the seminal papers [10–12] directed toward generalizing the dynamics of signal decay due to diffusion and due to relaxation (T1 and T2 ). In the case of diffusion, the connection is a simple extension of ordinary diffusion – with its linear relationship between the mean squared particle displacement and diffusion time (via D, the diffusion constant) – to anomalous diffusion where the diffusion time is raised to the fractional power (alpha, typically between zero and one). In relaxation the link is through generalization of the single exponential relaxation of homogeneous materials to the distribution of relaxation times, usually observed in complex materials, such as biological tissues. The fractional order model for NMR relaxation – as is the case for dielectrics and viscoelastic materials – is embedded in the Mittag-Leffler function, which for small arguments is well approximated by the stretched exponential. Both functions are easily shown to represent a distribution of exponential relaxation times [24]. The form of the solutions obtained in this study, of course, converged to the results of this previous work when the delay (tau) is zero as well as the classical solutions to NMR precession and relaxation when the fractional order of the derivative is set to equal one. One aspect of the fractional order model is the interconnection between precession and relaxation that arises due to the fact that, unlike the ordinary exponential, the Mittag-Leffler function cannot be simply factored. In fact, this model is very much like that reported in [25,26] where solutions to the damped fractional oscillator are expressed by them in terms of fractional exponential functions via fractional trigonometry. In the future, experimental NMR work is needed to search for such behavior. Additionally, other fractional order generalizations of the Bloch equation should be explored that do not affect precession. The application of delay to the fractional order Bloch equation is new, and to our knowledge, unique to this paper. In fact, although delays are a common currency of NMR experiments, NMR system models with intrinsic delay have not been encountered previously by us. Differential equations with delay, however, have almost as long a history as that of fractional calculus — dating back at least to Laplace and Condorcet in the eighteenth century [27]. Today, differential equations with delay are a mainstay of population biology through the so-called predator–prey equations that predict future population growth based not on current numbers of a given species, but on the population sizes at some time in the past (see, for example Refs. [28,29]). The solutions to such problems are not easy, even in the simplest cases, as one must specify an initial history function and test for the stability of the solutions [27]. Here we have followed the recent techniques from [27] to solve the fractional order model of the Bloch equation with delay. In our analysis, we distinguished two subsystems (A for T2 relaxation, and B for T1 relaxation). In the case of T1 relaxation, our Eq. (4.5) is equivalent to Eq. (10) analyzed in [30,31] and Eq. (7) analyzed in [17]. These two groups found that this subsystem is stable when Kp is less than one. Since, in our model, Kp is minus one, we anticipate stable behavior and that in fact was what was observed in our numerical examples for delays on the order of 10–100 µs. For much shorted delays (order of 1 µs) both subsystems A and B were stable, and mimicked typical fractional order relaxation expressed through the Mittag-Leffler function (with a complex argument). As the delay increased for subsystem A (T2 relaxation), this initial relaxation response was lost and the system transitioned to an unstable state with growing oscillations of the transverse magnetization. The exact value depended upon the fractional order. It is, of course, also dependent upon the precessional frequency (160 Hz) and the T2 (20 ms). Applying the stability analysis obtained in [17] gives in all cases a reasonable estimate of this system transition. One other difference between the two systems was the initial history function, g (t ), that is needed to solve the delay equations. For the Mz (t ) (T1 relaxation case), g (t ) must be a non-zero constant (we assumed 100) for t < 0, while for Mx (t ) and My (t ) (T2 relaxation case), g (t ) must be zero for t < 0. Immediately after time zero when we have assumed that a perfect (π /2)x pulse was applied, we have the initial conditions: Mx (0) = 0, My (0) = 100, and Mz (0) = 0. After a long time (many T1 s), the x and y magnetization components relax to zero while the z component returns to its equilibrium values — here assumed to be 100. 7. Conclusion The results presented in this paper demonstrate that solutions to the generalized fractional order Bloch equation with delay exhibit novel behavior dependent upon the selected system parameters. In general, for short delays the system relaxes with fractional order – not exponential – decay. While the T1 relaxation is stable for all values of the chosen delay (τ ), T2 relaxation is not and after a critical value of τ exhibits an unstable and increasing sinusoid oscillation. The emergence of this behavior is consistent with the current stability theory for these systems. The likelihood of such behavior occurring in experimental systems is unknown at present, but possible, as such growing oscillations (spin echoes) are observed in integer order systems following multiple RF excitation pulses. Further work is needed to connect the new fractional order models with delay to the NMR behavior of complex molecules and materials. Acknowledgements R.L. Magin would like to acknowledge the support of NIH grant R01 EB 007537 from NIBIB. V. Daftardar-Gejji acknowledges University of Pune, India and Department of Science and Technology, N. Delhi, India for the Research Grants. References [1] A. Abragam, Principles of Nuclear Magnetism, Oxford University Press, New York, 2002. [2] E.M. Haacke, R.W. Brown, M.R. Thompson, R. Venkatesan, Magnetic Resonance Imaging: Physical Principles and Sequence Design, John Wiley and Sons, New York, 1999. [3] Z.P. Liang, P.C. Lauterbur, Principles of Magnetic Resonance Imaging: A Signal Processing Perspective, IEEE Press, New York, 2000.
S. Bhalekar et al. / Computers and Mathematics with Applications 61 (2011) 1355–1365
1365
[4] M.A. Bernstein, K.F. King, X.J. Zhou, Handbook of MRI Pulse Sequences, Elsevier Academic Press, Burlington, MA, 2004. [5] M.T. Vlaardingerbroek, J.A. den Boer, Magnetic Resonance Imaging, second ed., Springer-Verlag, Berlin, 1999. [6] J. Kärger, G. Vojta, On the use of NMR pulsed field-gradient spectroscopy for the study of anomalous diffusion in fractal networks, Chem. Phys. Lett. 141 (1987) 411–413. [7] J. Kärger, H. Pfeifer, G. Vojta, Time correlation during anomalous diffusion in fractal systems and signal attenuation in NMR field-gradient spectroscopy, Phys. Rev. A 37 (1988) 4514–4517. [8] A. Widom, H.J. Chen, Fractal Brownian motion and nuclear spin echoes, J. Phys. A 28 (1998) 1243–1247. [9] R. Kimmich, Strange kinetics, porous media, and NMR, Chem. Phys. 284 (2002) 253–285. [10] A.E. Sitnitsky, G.G. Pimenov, A.V. Anisimov, Spin–lattice NMR relaxation by anomalous translational diffusion, J. Magn. Reson. 172 (1) (2005) 48–55. [11] R.L Magin, O. Abdullah, D. Baleanu, X.H.J. Zhou, Anomalous diffusion expressed through fractional order differential operators in the Bloch–Torrey equation, J. Magn. Reson. 190 (2) (2008) 255–270. [12] I. Petras, Modeling and numerical analysis of fractional-order Bloch equations, Comput. Math. Appl. 61 (2) (2010) 341–356. [13] R.L. Magin, X. Feng, D. Baleanu, Solving the Fractional Order Bloch Equation, Concept. Magn. Reson. Part A 34A (1) (2009) 16–23. [14] A.A. Kilbas, H.M. Srivastava, J.J. Trujillo, Theory and Applications of Fractional Differential Equations, in: North-Holland Mathematics Studies, vol. 204, 2006. [15] I. Podlubny, Fractional Differential Equations, Academic Press, San Diego, 1999. [16] S.G. Samko, A.A. Kilbas, O.I. Marichev, Fractional Integrals and Derivatives: Theory and Applications, Gordon and Breach, Yverdon, 1993. [17] W. Deng, C. Li, J. Lu, Stability analysis of linear fractional differential system with multiple time delays, Nonlinear Dynam. 48 (2007) 409–416. [18] K. Diethelm, N.J. Ford, A.D. Freed, A predictor–corrector approach for the numerical solution of fractional differential equations, Nonlinear Dynam. 29 (2002) 3–22. [19] K. Diethelm, An algorithm for the numerical solution of differential equations of fractional order, Elecron. Trans. Numer. Anal. 5 (1997) 1–6. [20] C.P. Li, C.X. Tao, On the fractional Adams method, Comput. Math. Appl. 58 (8) (2009) 1573–1588. [21] S. Bhalekar, V. Daftardar-Gejji, Chaos in nonlinear delay differential equations of fractional order, Comput. Math. Appl. (submitted for publication). [22] I. Podlubny, Matrix approach to discrete fractional calculus, Fract. Calc. Appl. Anal. 3 (4) (2000) 359–386. [23] I. Podlubny, A. Chechkin, T. Skovranek, Y.Q. Chen, B.M. Vinagre Jara, Matrix approach to discrete fractional calculus II: partial fractional differential equations, J. Comput. Physics (2009). [24] M.N. Berberan-Santos, B. Valeur, Luminescence decays with underlying distributions: general properties and analysis with mathematical functions, J. Lumin. 126 (2007) 263–272. [25] C.F. Lorenzo, T.T. Hartley, Fractional trigonometry and the spiral functions, Nonlinear Dynam. 38 (1-4) (2004) 23–60. [26] C.F. Lorenzo, T.T. Hartley, The fractional hyperbolic functions: with application to fractional differential equations, in: Proceedings of the ASME International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, vol. 6, Pts A-C, 2005, pp. 1505–1512. [27] H. Gorecki, S. Fuksa, P. Gabrowski, A. Koritowski, Analysis and Synthesis of Time Delay Systems, J. Willey and Sons-PWN, Warsaw, 1989. [28] F.C. Hoppensteadt, Z. Jackiewicz, Numerical solution of a problem in the theory of epidemics, Appl. Numer. Math. 56 (2006) 533–543. [29] A. Martin, S. Ruan, Predator–prey models with delay and prey harvesting, J. Math. Biol. 43 (2001) 247–267. [30] Y.Q. Chen, K.L. Moore, Analytical stability bound for a class of delayed fractional-order dynamic systems, Nonlinear Dynam. 29 (2002) 191–200. [31] R.D. Driver, Ordinary and Delay Differential Equations, in: Appl. Math. Sci., vol. 20, Springer-Verlag, New-York, 1977.
Contents lists available at ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
Fractional Bloch equation with delay Sachin Bhalekar a , Varsha Daftardar-Gejji a , Dumitru Baleanu b,c , Richard Magin d,∗ a
Department of Mathematics, University of Pune, Pune - 411007, India
b
Department of Mathematics and Computer Science, Faculty of Arts and Sciences, Cankaya University - 06530, Ankara, Turkey
c
Institute of Space Sciences, P.O. Box, MG-23, R 76900, Magurele-Bucharest, Romania
d
Department of Bioengineering, University of Illinois, 851 S. Morgan St., Chicago, 60607, USA
article
abstract
info
Article history: Received 10 July 2010 Received in revised form 30 December 2010 Accepted 31 December 2010
In this paper we investigate a fractional generalization of the Bloch equation that includes both fractional derivatives and time delays. The appearance of the fractional derivative on the left side of the Bloch equation encodes a degree of system memory in the dynamic model for magnetization. The introduction of a time delay on the right side of the equation balances the equation by also adding a degree of system memory on the right side of the equation. The analysis of this system shows different stability behavior for the T1 and the T2 relaxation processes. The T1 decay is stable for the range of delays tested (1–100 µs), while the T2 relaxation in this model exhibited a critical delay (typically 6 µs) above which the system was unstable. Delays are expected to appear in NMR systems, in both the system model and in the signal excitation and detection processes. Therefore, by including both the fractional derivative and finite time delays in the Bloch equation, we believe that we have established a more complete and more realistic model for NMR resonance and relaxation. © 2011 Elsevier Ltd. All rights reserved.
Keywords: Fractional calculus Bloch equation Delay
1. Introduction The classical description of nuclear magnetic resonance (NMR) – the phenomena underlying magnetic resonance imaging (MRI) – is summarized in vector form by the Bloch equation [1,2], →
dM dt
→
→
=γ M× B −
(Mz − M0 ) T1
iz −
(Mxix + Myiy ) T2
,
→
(1.1) →
where M (Mx , My , Mz ) represents the time-varying magnetization (Amps/s), M0 the equilibrium magnetization, B (Bx , By , Bz ) the applied radio frequency (RF), gradient and static magnetic fields (Tesla), γ is the gyromagnetic ratio (42.57 MHz/T, for spin 1/2 protons), and T1 , T2 are the spin–lattice and the spin–spin relaxation times, respectively. For homogeneous (e.g., over a 1–2 mm3 voxel), and isotropic materials with a single spin component (typically water protons), the Bloch equation prescribes the dynamic balance between externally applied magnetic fields and internal sample relaxation times. This dynamic balance is the basis for pulse sequence design, signal acquisition, image reconstruction and, in the case of MRI, tissue contrast [3–5]. Fractional order generalization of the Bloch equation has been undertaken by several groups [6–13] to account for the anomalous relaxation and anomalous diffusion observed in NMR studies of complex materials — typically gels, emulsions,
∗
Corresponding author. Tel.: +1 312 413 5528; fax: +1 312 996 5921. E-mail addresses: [email protected] (S. Bhalekar), [email protected], [email protected] (V. Daftardar-Gejji), [email protected] (D. Baleanu), [email protected] (R. Magin). 0898-1221/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.camwa.2010.12.079
1356
S. Bhalekar et al. / Computers and Mathematics with Applications 61 (2011) 1355–1365
porous composites and biological tissues. A common feature of such complex materials is a ‘‘mesoscopic’’ structure intermediate in scale between the molecular and the macroscopic regimes. NMR spectroscopy and MRI are powerful experimental tools for probing the organization (ordered/disordered) and the dynamics (phase transitions, diffusion, permeability) of mesoscopic structures. High resolution MR microscopy, for example, provides high contrast images of phase separation, particle aggregation, and macromolecular coalescence, while NMR diffusion studies give a direct measure of material porosity and tortuosity. However, the application of these experimental tools and the analysis of the acquired data are highly dependent upon the theoretical assumptions underlying the Bloch equation: typically Gaussian spin dynamics, Brownian particle motion, and first order exponential relaxation. To the extent that these assumptions apply to a particular material or to a transition between different material phases, conventional NMR analysis is valid and appropriate. Nevertheless, accumulating experimental evidence on complex materials, acquired using dielectric spectroscopy, dynamic light scattering as well as ultrasonic and viscoelastic measurements, strongly suggests anomalous dynamic behavior. This anomalous behavior appears to reflect distributions of relaxation times, and multi-scale phenomena that in some cases suggest a fractallike structure, nonlocal interactions, fading memory and aging. Hence, it is anticipated that NMR and MRI measurements of complex materials with mesoscopic states or configurations will show similar – fractional order – dynamic behavior. The largely unresolved question concerning the fractional or ‘‘generalized’’ Bloch equation is how to modify it to relate the fractional order of the operators to the spin dynamics at the molecular, mesoscopic, and macroscopic scales. NMR is a widely used experimental technique because it effectively probes time scales from the nanosecond (spin precession) to the microsecond (molecular reorientation) to the millisecond (spin–spin dephasing) and up to the second (spin–lattice relaxation) and beyond (diffusion). Using images or localized spectroscopy, NMR analysis assigns to each region of interest (pixel or voxel) specific experimental values (e.g., diffusion coefficient, spin–spin and spin–lattice relaxation times, correlation coefficient, chemical shift). This process works well for imaging macroscopic regions (i.e., voxels on the order of 1 mm3 ) or for spectroscopy on pure samples in uniform magnetic fields (i.e., one and two dimensional NMR spectra). However, pushing the imaging resolution into the microscopic or mesoscopic regimes or pulling complex spectra from macroscopic mixtures or tissues exposes fundamental limitations in either signal-to-noise (SNR) or spectral resolution. These problems, unfortunately, cannot always be overcome by taking more signal averages or by the use of extremely high field NMR systems. The goal of fractional order NMR modelling is to extend the experimental window of MR techniques in both space and time, not by simply increasing system resolution or bandwidth, but by appending to the governing equations a fractional order interpretation for the observed phenomena. The power of fractional calculus is that it conveniently assimilates into the order of the operator a fractal network of lacunae extending over a wide distribution of relaxation times or diffusion coefficients. Such conciseness comes at the cost of relinquishing the specific sub-compartment values for individual relaxation times or diffusion constants, but provides the benefit of capturing hierarchical, multi-scale phenomena in space and time. In previous work we have investigated spatial complexity via fractional order space derivatives [11] and anomalous temporal relaxation through fractional order time derivatives [12]. Both approaches appear promising for evaluating porous, layered or semi-permeable composite materials. In this paper we investigate – for the case of relaxation – the effect of introducing a short time delay into the fractional order model. Such time delays may perhaps be used to capture some of the longer time dynamic changes in complex systems associated with aging in both the experimental and system senses. In this paper we investigate – for the case of relaxation – the effect of introducing a short time delay into the fractional order model. Such time delays are common in NMR pulse programming where they are used to generate coherence in the detected signals, and to account for imperfections in the applied gradients and radio frequency pulses. Such time delays may also capture some of the longer time dynamic changes in complex systems associated with aging in both the experimental and system senses. Fractional order operators in time are well known to incorporate system memory; this arises through the inclusion of a power law function in the time domain convolution operator definition of fractional order integration and differentiation. Fractional order operators in space – in a similar manner – extend spatial derivatives through a non-local or distributed spatial gradient that captures a wider range of sample structure in its formal definition. The delay parameter investigated in this paper modifies the time window for the fractional order time derivative by allowing not only fading memory of an earlier state or time, but by also introducing specific information about the initial conditions at a time specified in the lower limit of the convolution operator. The rationale for this approach is that we hope the generalization will extend the modelling capabilities of the Bloch equation so that transient mesoscopic structures in complex materials can be better visualized. The organization of this paper is as follows: in Section 2 (Preliminaries) we present the fundamental definitions for the fractional order operators to be used; in Section 3 (Methods) we outline the numerical method used to solve the fractional Bloch equation with delay; in Section 4 (Theory) we describe the fractional Bloch equation with both fading memory and delay included; in Section 5 (Results) we present our findings for examples selected to span a range of fractional orders and short time delay; and in the last two sections we summarize our contribution by comparing and contrasting our results with our previous work and the work of others. 2. Preliminaries Basic definitions and properties of fractional derivative/integrals are given below [14–16].
S. Bhalekar et al. / Computers and Mathematics with Applications 61 (2011) 1355–1365
1357
Definition 2.1. A real function f (t ), t > 0 is said to be in space Cα , α ∈ ℜ if there exists a real number p (>α), such that f (t ) = t p f1 (t ) where f1 (t ) ∈ C [0, ∞). Definition 2.2. A real function f (t ), t > 0 is said to be in space Cαm , m ∈ N
{0} if f (m) ∈ Cα .
Definition 2.3. Let f ∈ Cα and α ≥ −1, then the (left-sided) Riemann–Liouville integral of order µ, µ > 0 is given by µ
I f (t ) =
t
∫
1
Γ (µ)
(t − τ )µ−1 f (τ )dτ ,
t > 0.
(2.1)
0
m Definition 2.4. The (left sided) Caputo fractional derivative of f , f ∈ C− 1, m ∈ N
Dµ f (t ) =
dm dt m
f (t ),
= I m−µ
{0}, is defined as:
µ=m
dm f ( t ) dt m
,
m − 1 < µ < m, m ∈ N .
(2.2)
Note that for m − 1 < µ ≤ m, m ∈ N, I µ Dµ f (t ) = f (t ) −
m−1
− dk f k =0
I µt ν =
tk
(0) , k!
dt k
(2.3)
Γ (ν + 1) t µ+ν . Γ (µ + ν + 1)
(2.4)
Stability analysis Deng et al. [17] have studied the stability of n-dimensional linear fractional differential equation with time delays for the following system: Dα1 x1 (t ) = a11 x1 (t − τ11 ) + a12 x2 (t − τ12 ) + · · · + a1n xn (t − τ1n )
Dα2 x2 (t ) = a21 x1 (t − τ21 ) + a22 x2 (t − τ22 ) + · · · + a2n xn (t − τ2n )
··· Dαn xn (t ) = an1 x1 (t − τn1 ) + an2 x2 (t − τn2 ) + · · · + ann xn (t − τnn ),
(2.5)
where αi ∈ (0, 1), xi (t ), xi (t − τij ) state variables and τij are delay terms. It can be shown that if all the roots of the characteristic equation det(∆(λ)) = 0, where
α λ 1 − a11 e−λτ11 −a21 e−λτ21 ∆(λ) = . ..
−a12 e−λτ12 λα2 − a22 e−λτ22 .. .
−an1 e−λτn1
−an2 e−λτn2
··· ··· .. . ···
−a1n e−λτ1n −a2n e−λτ2n .. .
λαn − ann e−λτnn
,
(2.6)
have negative real parts, then the system (2.5) is asymptotically stable. 3. Algorithm for solving delay differential equations of fractional order Due to the nonlocal nature of fractional derivative operator, the numerical methods for solving ordinary differential equations have to be modified to cater for solving fractional differential equations (FDEs) [18–20]. Bhalekar and DaftardarGejji have modified Adams–Bashforth method to solve delay differential equations of fractional order (FDDE) [21]. The method is described below. Consider the FDDE Dα y(t ) = f (t , y(t − τ )), y(t ) = g (t ),
t ∈ [0, T ], 0 < α ≤ 1,
t ∈ [−τ , 0],
(3.1) (3.2)
where Dα denotes the Caputo fractional derivative of order α ∈ (0, 1]. Applying fractional integration Itα on both sides of (3.1) and using (3.2) we obtain y(t ) = g (0) +
1
Γ (α)
t
∫ 0
(t − ξ )α−1 f (ξ , y(ξ − τ ))dξ .
(3.3)
1358
S. Bhalekar et al. / Computers and Mathematics with Applications 61 (2011) 1355–1365
Consider a uniform grid {tn = nh : n = −k, −k + 1, . . . , −1, 0, 1, . . . , N } where k and N are integers such that h = T /N and h = τ /k. Let yh (tj ) = g (tj ),
j = −k, −k + 1, . . . , −1, 0
(3.4)
and note that yh (tj − τ ) = yh (jh − kh) = yh (tj−k ),
j = 0, 1, . . . , N .
(3.5)
Suppose we have already calculated approximations yh (tj ) ≈ y(tj ), (j = −k, −k + 1, . . . , −1, 0, 1, . . . , n) and we want to calculate yh (tn+1 ) using tn+1
∫
1
(tn+1 − ξ )α−1 f (ξ , y (ξ − τ ))dξ . (3.6) Γ (α) 0 We use approximations yh (tn ) for y(tn ) in (3.6). The integral on right hand side of (3.6) is evaluated using the product y(tn+1 ) = g (0) +
trapezoidal quadrature formula. We obtain yh (tn+1 ) = g (0) +
= g (0) +
hα
Γ (α + 2) hα
Γ (α + 2)
f (tn+1 , yh (tn+1 − τ )) +
f (tn+1 , yh (tn+1−k )) +
hα
n −
Γ (α + 2)
j =0
hα
n −
Γ (α + 2)
j =0
aj,n+1 f (tj , yh (tj − τ ))
aj,n+1 f (tj , yh (tj−k )),
(3.7)
where aj,n+1 are given by aj,n+1
nα+1 − (n − α)(n + 1)α , = (n − j + 2)α+1 + (n − j)α+1 − 2(n − j + 1)α+1 , 1,
if j = 0, if 1 ≤ j ≤ n, if j = n + 1.
(3.8)
Note that, unlike in the non-delayed case the right hand side of the Eq. (3.7) does not contain the unknown term yh (tn+1 ), so it is not needed to calculate the predictor term. 4. Fractional Bloch equation with delay We apply the method described in Section 3 for solving fractional order Bloch equation involving time delay T −α Dα Mx (t ) = ω˜0 My (t − τ ) −
Mx (t − τ ) T2
T −α Dα My (t ) = −ω˜0 Mx (t − τ ) − T −α Dα Mz (t ) = Mx (t ) = 0, where ω˜0 =
ω0
T α−1
M0 − Mz (t − τ ) T1 My (t ) = 100,
= (ω0 T )T −α ,
1 T′
=
1
,
My (t − τ ) T2
(4.1)
,
(4.2)
,
(4.3) Mz (t ) = 0,
1 T α−1 T
= 1
for t ≤ 0,
T −α 1 T , T′ T1 2
=
(4.4)
1 T α−1 T
2
=
T −α T . We assume that T T2 1 1 0 T ′ and T ′ have the unit
ω˜0 = 3 × π , T1′ = 1, T2′ = 20 × 10−3 , M0 = 100. We mention that ω˜ ,
1
2
= 1.0. In our case we have of (s)−α .
We use the stability analysis given in Section 2 to study the behavior of the system (4.1)–(4.3). Note that, Eq. (4.3) can be transformed to Dα Mz1 (t ) = −Mz1 (t − τ ),
(4.5)
where Mz1 (t ) = (Mz (t ) − M0 )/T1 . The characteristic equation for this system becomes:
λ2α +
2 α −λτ λ e + T2
1
1
e−2λτ
(λα + e−λτ ) = 0. (4.6) ω02 T22 + 100λα e−λτ + (102400π 2 + 2500)e−2λτ ) in (4.6) corresponds to Eq. (4.1)–(4.2) and the factor +
Note that the factor (λ2α (λα + e−λτ ) corresponds to a decoupled Eq. (4.5). We rename (4.1)–(4.2) as subsystem (A) and (4.5) as subsystem (B). The system changes its behavior (from damped oscillations to oscillations with an increasing height) for a suitable value of delay. Note that when the delay is zero, the solutions to this system of fractional order differential equations are given in [12]. When α = 1 and the delay is zero, we find the classical results. 5. Results The observations are summarized in Figs. 1–6. We consider three values of α for both system A and system B. In this example we assume a perfect spin excitation ( π2 )x so that Mx (0) = 0, My (0) = M0 , Mz (0) = 0. This system should return to equilibrium at long times compared with T1 and T2 , that is Mx (∞) = 0, My (∞) = 0, Mz (∞) = M0 .
S. Bhalekar et al. / Computers and Mathematics with Applications 61 (2011) 1355–1365
a
b
Mx
Mz
c
My 100
80
50
60
100 50 t 0.01
0.02
0.03
0.01
0.04
0.03
0.04
40
t
20
-50
-50
-100
-100
t 1
My
d
0.02
1359
2
3
4
e
100
50
-100
-50
50
100
Mx
100 50
M z 642 0 -100
-50
0 M y -50
-50
0 -100
Mx
50
-100
100
Fig. 1. A set of plots of magnetization versus time for the Mx , My and Mz relaxation dynamics following a 90x degree pulse excitation (a, b, and c respectively) and a two-dimensional plot of Mx versus My (d) and a three-dimensional plot of Mx , My and Mz (e). In this figure the time delay was 0.00006 s and the value of α was 0.8.
a
b
Mx
150000
150000
100000
100000
50000
50000 0.01
0.02
0.03
0.04
c
-100000
Mz 80
t
-50000
0.01
0.02
0.03
0.04
t
60
-50000
40
-100000
20
-150000
-150000
d
My
t 1
My
2
3
4
e
150000 100000 50000 Mx
-150000 -100000 -50000
50000 100000 150000 2000
-50000 -100000 -150000
Mz 6402 -20000
0 My 0 Mx
-20000 20000
Fig. 2. A set of plots of magnetization versus time for the Mx , My and Mz relaxation dynamics following a 90x degree pulse excitation (a, b, and c respectively) and a two-dimensional plot of Mx versus My (d) and a three-dimensional plot of Mx , My and Mz (e). In this figure the time delay was 0.00007 s and the value of α was 0.8.
In the case of α = 0.80: for τ = 0.00006, roots of the characteristic Eq. (4.6) corresponding to system (A) are −143.972 ± 5729.25i and corresponding to system (B) are −163062 ± 123853i. In view of the analysis given in Section 2,
1360
S. Bhalekar et al. / Computers and Mathematics with Applications 61 (2011) 1355–1365
a
c
b
Mx
100
100
80
50
60
50 0.01
0.02
0.03
t 0.04
0.01
-50
-50
-100
-100
My
d
Mz
My
0.02
40
t 0.04
0.03
20 t 1
2
3
4
e
100
50
-100
-50
50
Mx
100
100 50
M z 43210
-50
0 M y
-100 -50
-50
0 -100
50
Mx
-100 100
Fig. 3. A set of plots of magnetization versus time for the Mx , My and Mz relaxation dynamics following a 90x degree pulse excitation (a, b, and c respectively) and a two-dimensional plot of Mx versus My (d) and a three-dimensional plot of Mx , My and Mz (e). In this figure the time delay was 0.00009 s and the value of α was 0.9.
a
b
Mx
c
My
200
200
100
100
Mz 80
0.01
0.02
0.03
0.04
60
t
0.01
-100
-100
-200
-200
0.02
0.03
0.04
t
40 20 t 1
2
3
4
My
d
e
200
100
-200
-100
100
-100
200
Mx 200 100
Mz 43210 -200
0 My
-100 Mx -200
-100
0 100 200
-200
Fig. 4. A set of plots of magnetization versus time for the Mx , My and Mz relaxation dynamics following a 90x degree pulse excitation (a, b, and c respectively) and a two-dimensional plot of Mx versus My (d) and a three-dimensional plot of Mx , My and Mz (e). In this figure the time delay was 0.00010 s and the value of α was 0.9.
the system has stable solutions in this case. The numerical results also have the same conclusions. The damped oscillations of Mx and My are shown in Fig. 1(a) and (b), whereas Fig. 1(c) shows the stable solution for Mz .
S. Bhalekar et al. / Computers and Mathematics with Applications 61 (2011) 1355–1365
1361
c a
b
Mx
Mz 100
My
100
100
80
50
50
60
0.01
0.02
0.03
t
0.04
0.01
-50
-50
-100
-100
0.02
0.03
0.04
40
t
20 1
2
3
t
4
My
d
e
100
50
-50
50
100
Mx 100 50
Mz 3210
-50
0 My
-50 Mx
-50
0 50
100 -100
-100
Fig. 5. A set of plots of magnetization versus time for the Mx , My and Mz relaxation dynamics following a 90x degree pulse excitation (a, b, and c respectively) and a two-dimensional plot of Mx versus My (d) and a three-dimensional plot of Mx , My and Mz (e). In this figure the time delay was 0.00004 s and the value of α was 1.0.
b
a
Mx
c
My 150
150
80
100
100
60
50
50 0.01
0.02
0.03
0.04
t
0.01
-50
-50
-100
-100
-150
-150
d
Mz 100
0.02
0.03
0.04
t
40 20 1
2
3
4
t
e
My 150 100 50
-150
-100
-50
50
100
150
Mx
-50
10 0 Mz 43210
-100 -150
0 My
-100 Mx
0
-100 100
Fig. 6. A set of plots of magnetization versus time for the Mx , My and Mz relaxation dynamics following a 90x degree pulse excitation (a, b, and c respectively) and a two-dimensional plot of Mx versus My (d) and a three-dimensional plot of Mx , My and Mz (e). In this figure the time delay was 0.00006 s and the value of α was 1.0.
1362
S. Bhalekar et al. / Computers and Mathematics with Applications 61 (2011) 1355–1365
When τ is increased to 0.00007, the roots of characteristic equation become 184.389 ± 5575.72i (system (A)) and −135039 ± 9804.27i (system (B)). Since the real parts of roots for system (A) are non-negative, we cannot use the theory in Section 2. The numerical results shows that the system (A) becomes unstable for these values of parameters. In Fig. 2(a) and (b) the oscillations of Mx and My are with increasing height. This behavior is not normally observed in NMR experiments without delay unless addition RF pulses are applied that generate so-called spin echoes. The system (B) is stable (Fig. 2(c)) in this case. In the case of α = 0.90: damped oscillations of Mx and My are observed for τ = 0.00009 (Fig. 3(a), (b)). This observation by numerical methods is supported by the theory also, as the roots of the characteristic equation corresponding to system (A) are −26.4959 ± 2175.72i. Also for system (B), the roots are −116677 ± 3817.75i and the system is stable (Fig. 3(c)). For τ = 0.00010, the characteristic roots become 23.2458 ± 2164.4i (system (A)) and −103972 ± 3439.19i (system (B)). The numerical solutions (Fig. 4(a),(b)) shows that the system (A) is unstable and system (B) (Fig. 4(c)) is stable. In the case of α = 1.00: it can be observed from analytically and numerically that the system is stable for τ = 0.00004 (Fig. 5(a)–(c)). This is the expected classical behavior. For τ = 0.00006, oscillations of Mx , My with increasing height are observed (Fig. 6(a), (b)) and Mz is stable (Fig. 5(d)). Thus, even for the integer order time derivative, a short time delay is sufficient to bring about a shift in the dynamical behavior of the system. In addition all three values of α(0.80, 0.90, 1.0) a very short time delay (τ = 0.000001 s) yielded only stable anomalous relaxations via Mittag-Leffler functions. In this example we have assumed ω0 = 160 × 2π rad/s. This corresponds to a sinusoidal oscillator with a period of (T0 = 2ωπ ) of 6.25 ms, while the delay typically ranges from 0 40–100 µs. No available behavior was observed for the delay of 1 µs. 5.1. Comparison with existing method In this subsection, we compare our results with the existing discretization method discussed in [15,22,23] by Podlubny and co-workers. Let f (t ) be a function defined on [0, b]. Consider the equidistant nodes ti = ih (i = 0, 1, . . . , N ) with the step h, where t0 = 0 and tN = b. If f (j) (0) = 0, j = 0, 1, . . . , m then we can approximate fractional order derivative Dα f (t ) of order α (m − 1 < α ≤ m) at point ti as follows: Dα f (t ) |t =ti ≈ h−α
i α − (−1)r fi−r ,
(5.1)
r
r =0
where
α r
=
Γ (α + 1) r !Γ (α − r + 1)
and fi = f (ti ). Using the transformation x(t ) = Mx (t ), y(t ) = My (t ) − 100, z (t ) = Mz (t ), we can write the system (4.1)–(4.3) as Dα x(t ) = ω˜0 (y(t − τ ) + 100) −
M0 T1
x(t ) = 0,
−
z (t − τ ) T1
T2
y(t − τ ) + 100
Dα y(t ) = −ω˜0 x(t − τ ) − Dα z (t ) =
x(t − τ )
T2
,
(5.2)
,
(5.3)
,
y(t ) = 0,
(5.4) z (t ) = 0,
for t ≤ 0, 0 < α ≤ 1.
(5.5)
Consider a uniform grid {tn = nh : n = −k, −k + 1, . . . , −1, 0, 1, . . . , N } where k and N are integers such that h = b/N and h = τ /k. Using (5.1) we can write the approximation of system (5.2)–(5.5) as h−α
i −
(−1)r
r =0
h−α
i −
(−1)r
r =0
h−α
i − r =0
(−1)r
α r
α r
α r
xi−r = ω˜0 (yi−k + 100) − yi−r = −ω˜0 xi−k −
zi−r =
M0 T1
−
zi−k T1
,
x i −k T2
yi−k + 100 T2
, ,
xi = yi = zi = 0, for i ≤ 0.
S. Bhalekar et al. / Computers and Mathematics with Applications 61 (2011) 1355–1365
a
b
Mx
1363
My
100 100 50 50
0.01
0.02
0.03
0.04
0.01
t
0.02
0.03
0.04
t
-50 -50 -100 -100
c
d
Mz
My 100
80 50 60
-100
40
-50
50
100
Mx
-50 20 -100 t 1
2
3
4
Fig. 7. Plots of magnetization versus time for the Mx , My and Mz relaxation dynamics following a 90x degree pulse excitation (a, b, and c respectively) and a two-dimensional plot of Mx versus My (d), with time delay 0.00006 s and α = 0.8 using the method given in Section 5.1.
Hence we can derive the following algorithm for obtaining a numerical solution xi = yi = zi = 0,
for i ≤ 0,
[
xi−k xi = hα ω˜0 (yi−k + 100) − T2
[
yi−k + 100 yi = hα −ω˜0 xi−k − T2 z i = hα
[
M0 T1
−
zi−k T1
] −
] −
i −
(−1)r
r =1
] −
r
xi−r ,
i α − (−1)r yi−r , r =1
i α − (−1)r zi−r , r =1
α
r
r
(5.6)
where i = 1, 2, . . . , N. In Fig. 7, we have plotted the case α = 0.8, τ = 0.0006 using the algorithm (5.6). Comparing these figures with Fig. 1, we can conclude that the results by our method agree well with the existing discretization method. From this we can conclude that our results match well with the classical method and are consistent with the stability analysis discussed in Section 2. 6. Discussion In this paper we have investigated a generalization of the Bloch equation that includes both fractional derivatives and time delays. The combination of these two powerful techniques is justified for by both theoretical and experimental reasons. First, in theory, the fractional derivative appears on the left side of the Bloch equation, hence the equation is asymmetric with the memory encoded in the definition of the fractional derivative contributing only to one part of the dynamic model for magnetization. In order to account for this discrepancy, we have introduced system memory on the right side of the Bloch equation through a finite time delay — tau. Second, experimentally, delays are expected to appear in NMR systems, delay in both the system model and in the signal excitation and detection processes. Therefore, by including both the fractional derivative and finite time delays in the Bloch equation, we believe that we have established a more complete and more realistic model for NMR resonance and relaxation. In the following we will describe more completely our rationale for this generalization and the implications of our results for the analysis of NMR signals. The application of fractional
1364
S. Bhalekar et al. / Computers and Mathematics with Applications 61 (2011) 1355–1365
calculus to NMR is relatively recent and it can be seen in the seminal papers [10–12] directed toward generalizing the dynamics of signal decay due to diffusion and due to relaxation (T1 and T2 ). In the case of diffusion, the connection is a simple extension of ordinary diffusion – with its linear relationship between the mean squared particle displacement and diffusion time (via D, the diffusion constant) – to anomalous diffusion where the diffusion time is raised to the fractional power (alpha, typically between zero and one). In relaxation the link is through generalization of the single exponential relaxation of homogeneous materials to the distribution of relaxation times, usually observed in complex materials, such as biological tissues. The fractional order model for NMR relaxation – as is the case for dielectrics and viscoelastic materials – is embedded in the Mittag-Leffler function, which for small arguments is well approximated by the stretched exponential. Both functions are easily shown to represent a distribution of exponential relaxation times [24]. The form of the solutions obtained in this study, of course, converged to the results of this previous work when the delay (tau) is zero as well as the classical solutions to NMR precession and relaxation when the fractional order of the derivative is set to equal one. One aspect of the fractional order model is the interconnection between precession and relaxation that arises due to the fact that, unlike the ordinary exponential, the Mittag-Leffler function cannot be simply factored. In fact, this model is very much like that reported in [25,26] where solutions to the damped fractional oscillator are expressed by them in terms of fractional exponential functions via fractional trigonometry. In the future, experimental NMR work is needed to search for such behavior. Additionally, other fractional order generalizations of the Bloch equation should be explored that do not affect precession. The application of delay to the fractional order Bloch equation is new, and to our knowledge, unique to this paper. In fact, although delays are a common currency of NMR experiments, NMR system models with intrinsic delay have not been encountered previously by us. Differential equations with delay, however, have almost as long a history as that of fractional calculus — dating back at least to Laplace and Condorcet in the eighteenth century [27]. Today, differential equations with delay are a mainstay of population biology through the so-called predator–prey equations that predict future population growth based not on current numbers of a given species, but on the population sizes at some time in the past (see, for example Refs. [28,29]). The solutions to such problems are not easy, even in the simplest cases, as one must specify an initial history function and test for the stability of the solutions [27]. Here we have followed the recent techniques from [27] to solve the fractional order model of the Bloch equation with delay. In our analysis, we distinguished two subsystems (A for T2 relaxation, and B for T1 relaxation). In the case of T1 relaxation, our Eq. (4.5) is equivalent to Eq. (10) analyzed in [30,31] and Eq. (7) analyzed in [17]. These two groups found that this subsystem is stable when Kp is less than one. Since, in our model, Kp is minus one, we anticipate stable behavior and that in fact was what was observed in our numerical examples for delays on the order of 10–100 µs. For much shorted delays (order of 1 µs) both subsystems A and B were stable, and mimicked typical fractional order relaxation expressed through the Mittag-Leffler function (with a complex argument). As the delay increased for subsystem A (T2 relaxation), this initial relaxation response was lost and the system transitioned to an unstable state with growing oscillations of the transverse magnetization. The exact value depended upon the fractional order. It is, of course, also dependent upon the precessional frequency (160 Hz) and the T2 (20 ms). Applying the stability analysis obtained in [17] gives in all cases a reasonable estimate of this system transition. One other difference between the two systems was the initial history function, g (t ), that is needed to solve the delay equations. For the Mz (t ) (T1 relaxation case), g (t ) must be a non-zero constant (we assumed 100) for t < 0, while for Mx (t ) and My (t ) (T2 relaxation case), g (t ) must be zero for t < 0. Immediately after time zero when we have assumed that a perfect (π /2)x pulse was applied, we have the initial conditions: Mx (0) = 0, My (0) = 100, and Mz (0) = 0. After a long time (many T1 s), the x and y magnetization components relax to zero while the z component returns to its equilibrium values — here assumed to be 100. 7. Conclusion The results presented in this paper demonstrate that solutions to the generalized fractional order Bloch equation with delay exhibit novel behavior dependent upon the selected system parameters. In general, for short delays the system relaxes with fractional order – not exponential – decay. While the T1 relaxation is stable for all values of the chosen delay (τ ), T2 relaxation is not and after a critical value of τ exhibits an unstable and increasing sinusoid oscillation. The emergence of this behavior is consistent with the current stability theory for these systems. The likelihood of such behavior occurring in experimental systems is unknown at present, but possible, as such growing oscillations (spin echoes) are observed in integer order systems following multiple RF excitation pulses. Further work is needed to connect the new fractional order models with delay to the NMR behavior of complex molecules and materials. Acknowledgements R.L. Magin would like to acknowledge the support of NIH grant R01 EB 007537 from NIBIB. V. Daftardar-Gejji acknowledges University of Pune, India and Department of Science and Technology, N. Delhi, India for the Research Grants. References [1] A. Abragam, Principles of Nuclear Magnetism, Oxford University Press, New York, 2002. [2] E.M. Haacke, R.W. Brown, M.R. Thompson, R. Venkatesan, Magnetic Resonance Imaging: Physical Principles and Sequence Design, John Wiley and Sons, New York, 1999. [3] Z.P. Liang, P.C. Lauterbur, Principles of Magnetic Resonance Imaging: A Signal Processing Perspective, IEEE Press, New York, 2000.
S. Bhalekar et al. / Computers and Mathematics with Applications 61 (2011) 1355–1365
1365
[4] M.A. Bernstein, K.F. King, X.J. Zhou, Handbook of MRI Pulse Sequences, Elsevier Academic Press, Burlington, MA, 2004. [5] M.T. Vlaardingerbroek, J.A. den Boer, Magnetic Resonance Imaging, second ed., Springer-Verlag, Berlin, 1999. [6] J. Kärger, G. Vojta, On the use of NMR pulsed field-gradient spectroscopy for the study of anomalous diffusion in fractal networks, Chem. Phys. Lett. 141 (1987) 411–413. [7] J. Kärger, H. Pfeifer, G. Vojta, Time correlation during anomalous diffusion in fractal systems and signal attenuation in NMR field-gradient spectroscopy, Phys. Rev. A 37 (1988) 4514–4517. [8] A. Widom, H.J. Chen, Fractal Brownian motion and nuclear spin echoes, J. Phys. A 28 (1998) 1243–1247. [9] R. Kimmich, Strange kinetics, porous media, and NMR, Chem. Phys. 284 (2002) 253–285. [10] A.E. Sitnitsky, G.G. Pimenov, A.V. Anisimov, Spin–lattice NMR relaxation by anomalous translational diffusion, J. Magn. Reson. 172 (1) (2005) 48–55. [11] R.L Magin, O. Abdullah, D. Baleanu, X.H.J. Zhou, Anomalous diffusion expressed through fractional order differential operators in the Bloch–Torrey equation, J. Magn. Reson. 190 (2) (2008) 255–270. [12] I. Petras, Modeling and numerical analysis of fractional-order Bloch equations, Comput. Math. Appl. 61 (2) (2010) 341–356. [13] R.L. Magin, X. Feng, D. Baleanu, Solving the Fractional Order Bloch Equation, Concept. Magn. Reson. Part A 34A (1) (2009) 16–23. [14] A.A. Kilbas, H.M. Srivastava, J.J. Trujillo, Theory and Applications of Fractional Differential Equations, in: North-Holland Mathematics Studies, vol. 204, 2006. [15] I. Podlubny, Fractional Differential Equations, Academic Press, San Diego, 1999. [16] S.G. Samko, A.A. Kilbas, O.I. Marichev, Fractional Integrals and Derivatives: Theory and Applications, Gordon and Breach, Yverdon, 1993. [17] W. Deng, C. Li, J. Lu, Stability analysis of linear fractional differential system with multiple time delays, Nonlinear Dynam. 48 (2007) 409–416. [18] K. Diethelm, N.J. Ford, A.D. Freed, A predictor–corrector approach for the numerical solution of fractional differential equations, Nonlinear Dynam. 29 (2002) 3–22. [19] K. Diethelm, An algorithm for the numerical solution of differential equations of fractional order, Elecron. Trans. Numer. Anal. 5 (1997) 1–6. [20] C.P. Li, C.X. Tao, On the fractional Adams method, Comput. Math. Appl. 58 (8) (2009) 1573–1588. [21] S. Bhalekar, V. Daftardar-Gejji, Chaos in nonlinear delay differential equations of fractional order, Comput. Math. Appl. (submitted for publication). [22] I. Podlubny, Matrix approach to discrete fractional calculus, Fract. Calc. Appl. Anal. 3 (4) (2000) 359–386. [23] I. Podlubny, A. Chechkin, T. Skovranek, Y.Q. Chen, B.M. Vinagre Jara, Matrix approach to discrete fractional calculus II: partial fractional differential equations, J. Comput. Physics (2009). [24] M.N. Berberan-Santos, B. Valeur, Luminescence decays with underlying distributions: general properties and analysis with mathematical functions, J. Lumin. 126 (2007) 263–272. [25] C.F. Lorenzo, T.T. Hartley, Fractional trigonometry and the spiral functions, Nonlinear Dynam. 38 (1-4) (2004) 23–60. [26] C.F. Lorenzo, T.T. Hartley, The fractional hyperbolic functions: with application to fractional differential equations, in: Proceedings of the ASME International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, vol. 6, Pts A-C, 2005, pp. 1505–1512. [27] H. Gorecki, S. Fuksa, P. Gabrowski, A. Koritowski, Analysis and Synthesis of Time Delay Systems, J. Willey and Sons-PWN, Warsaw, 1989. [28] F.C. Hoppensteadt, Z. Jackiewicz, Numerical solution of a problem in the theory of epidemics, Appl. Numer. Math. 56 (2006) 533–543. [29] A. Martin, S. Ruan, Predator–prey models with delay and prey harvesting, J. Math. Biol. 43 (2001) 247–267. [30] Y.Q. Chen, K.L. Moore, Analytical stability bound for a class of delayed fractional-order dynamic systems, Nonlinear Dynam. 29 (2002) 191–200. [31] R.D. Driver, Ordinary and Delay Differential Equations, in: Appl. Math. Sci., vol. 20, Springer-Verlag, New-York, 1977.