Physica Medica 32 (2016) 94–103

Contents lists available at ScienceDirect

Physica Medica j o u r n a l h o m e p a g e : h t t p : / / w w w. p h y s i c a m e d i c a . c o m

Original Paper

Improving the time efficiency of the Fourier synthesis method for slice selection in magnetic resonance imaging B. Tahayori a,*, N. Khaneja b, L.A. Johnston a, P.M. Farrell a, I.M.Y. Mareels a a b

Department of Electrical and Electronic Engineering, The University of Melbourne, Parkville, Victoria 3010, Australia School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA

A R T I C L E

I N F O

Article history: Received 11 July 2015 Received in revised form 6 October 2015 Accepted 9 October 2015 Available online 24 October 2015 Keywords: The Bloch equation Fourier synthesis Magnetic resonance imaging Slice selection Mathematical induction

A B S T R A C T

The design of slice selective pulses for magnetic resonance imaging can be cast as an optimal control problem. The Fourier synthesis method is an existing approach to solve these optimal control problems. In this method the gradient field as well as the excitation field are switched rapidly and their amplitudes are calculated based on a Fourier series expansion. Here, we provide a novel insight into the Fourier synthesis method via representing the Bloch equation in spherical coordinates. Based on the spherical Bloch equation, we propose an alternative sequence of pulses that can be used for slice selection which is more time efficient compared to the original method. Simulation results demonstrate that while the performance of both methods is approximately the same, the required time for the proposed sequence of pulses is half of the original sequence of pulses. Furthermore, the slice selectivity of both sequences of pulses changes with radio frequency field inhomogeneities in a similar way. We also introduce a measure, referred to as gradient complexity, to compare the performance of both sequences of pulses. This measure indicates that for a desired level of uniformity in the excited slice, the gradient complexity for the proposed sequence of pulses is less than the original sequence. © 2015 Associazione Italiana di Fisica Medica. Published by Elsevier Ltd. All rights reserved.

Introduction Magnetic resonance imaging (MRI) has become one of the most important medical imaging diagnostic tools available to physicians [1]. The basis of MRI is a phenomenon known as nuclear magnetic resonance (NMR). It is related to the way in which elementary particles such as protons interact with external magnetic fields having static and oscillating components. In MRI, it is possible to selectively excite the spins in a thin slice of the object by applying gradient fields [2]. To select a slice in a desired direction, a gradient field must be applied in that direction and a Radio Frequency (RF) excitation pulse with a limited bandwidth should be applied perpendicular to the gradient direction [2]. Many methods have been proposed to solve the selective excitation problem for adiabatic and non-adiabatic passages [3–9]. These methods are based on approximate solutions to the Bloch equation governing the spin dynamics in MRI, referred to as the Bloch equation [10–12], or are generated from computer simulations that predict the bulk magnetisation response to different excitation patterns [13,14]. Fourier analysis of the Bloch equation can be used to

* Corresponding author. Department of Electrical and Electronic Engineering, The University of Melbourne, Parkville, Victoria 3010, Australia. Tel.: +61 3 9035 3918; fax: +61 3 9035 3001. E-mail address: [email protected] (B. Tahayori).

find a selective pulse rotating the bulk magnetisation less than π/2, see [10] and [15]. In this method, it is assumed that the Bloch equation behaves linearly for small tip angles. The design of better pulses requires application of optimal control theory [16–19]. In [16] a mathematical basis for RF pulse design, an efficient algorithm to find the optimal pulse is provided. An optimal pulse is defined as the pulse that steers the magnetisation from the initial state, closest to the desired final state in a fixed amount of time. Ensemble controllability of the Bloch equation has been studied in [20,21]. The Shinnar–Le Roux (SLR) approach, which is ubiquitously used in designing pulses for MRI machines, is a recursive algorithm for finding the optimal pulse for a given selective excitation pattern [17]. In this approach the applied pulse is discretised to several rectangular pulses and the effect on the bulk magnetisation is calculated analytically at each step. As a result, the problem of finding selective pulses is reduced to the design of two polynomials. In this case, a selective RF pulse can be calculated through solving finite impulse response (FIR) filters. A number of variants of the SLR technique for designing selective pulses are available in the MRI literature [22–25]. Moore et al. in [5,7] have revisited the slice selection problem in the presence of the RF and the static field inhomogeneities using composite pulses. Their approach has two major steps. First, they optimise a non-selective composite pulse to minimise the effect of field inhomogeneities. In the second step, they use Gaussian or sinc sub-pulses with time varying gradient fields to enhance the slice

http://dx.doi.org/10.1016/j.ejmp.2015.10.088 1120-1797/© 2015 Associazione Italiana di Fisica Medica. Published by Elsevier Ltd. All rights reserved.

B. Tahayori et al./Physica Medica 32 (2016) 94–103

selectivity of the overall sequence of pulses. The uniformity of the bulk magnetisation flip angle has been improved in a 7 Tesla human scanner in the presence of field inhomogeneities. The Fourier synthesis method is an alternative technique to solve the optimal control problem for slice selection and RF field inhomogeneity suppression [26–30]. The gradient field as well as the excitation field are switched rapidly to shift the position of an ensemble of spins gradually. The magnitudes and duration of the applied fields are calculated based on a Fourier series expansion of the Bloch equation. In this paper, we propose a novel sequence of pulses for the Fourier synthesis method that considerably improves the time efficiency of this technique while preserving its efficiency in terms of the selected slice quality. Mathematical induction is employed to form a proof of the Fourier synthesis method in the spherical coordinates. Numerical simulation for both sequences of pulses is compared in this paper. Simulation results indicate that the slice selectivity of both sequences of pulses behaves similarly in the presence of RF field inhomogeneities.

The Bloch equation

1 ⎡ ⎤ Δω0 (r , t ) v (t ) ⎥ ⎢ − T (r)  2 ⎡M x ′ ⎤ ⎢ ⎥ ⎡M x ′ ⎤ ⎡ 0 ⎤ ⎢  ⎥ ⎢ ⎥ 1 ⎢ 1 u (t ) ⎥ ⎢M y ′ ⎥ + 0 ⎥, ⎢M y ′ ⎥ = ⎢ −Δω0 (r , t ) − ⎥ T1(r ) ⎢ ⎥ ⎢ T 2 (r) ⎢  ⎥ ⎢ ⎥ ⎢M ′ ⎥ ⎢⎣M 0 ⎥⎦ z ⎦ ⎣ ⎢⎣ M z ′ ⎥⎦ ⎢ 1 ⎥ −u (t ) − ⎢ −v (t ) ⎥ T1(r) ⎦ ⎣

(1)

u (t ) ≡ ωx ′ (t ) = γ B1e (t )cos φe ,

(2a)

v (t ) ≡ ωy ′ (t ) = γ B1e (t )sin φe ,

(2b)

and (3)

The parameters of the Bloch equation are summarised in Table 1. r is a vector representing position. If the excitation field is initially applied in the x-direction then ϕe = 0, and we may write

u (t ) ≡ ωx ′ (t ) ≡ ω1(t ) = γ B1e (t ), v (t ) ≡ ωy ′ (t ) = 0.

Parameter

Description

M(r,t)

Magnetisation vector in laboratory frame of reference Magnetisation vector in rotating frame of reference Bulk magnetisation magnitude at thermal equilibrium Longitudinal relaxation time constant Transverse relaxation time constant Gyromagnetic ratio External static field applied in the z-direction Excitation or rotating field applied in the xyplane Envelope of the excitation field Initial phase of the rotating field Gradient field Larmor frequency of the static field Oscillating frequency of the rotating field Rabi frequency Off-resonance excitation Deviation from the Larmor frequency as a result of the static field imperfection and the tiny fields induced by the object under study Space-dependent frequency generated by gradient fields Deviation from the Larmor frequency caused by all possible sources over space and time

M′(r,t) M0 T1(r) T2(r) γ B0 B1 B1e φe G r (t) ω0 = γB0 ωrf ω1(t) Δωoff δω0 = γδB0

Δω0 (r , t ) = Δωoff + δω0 + Δω (r , t )

where u(t) and v(t) are defined in Eq. (2). The above equation may be written as

⎡M x ′ ⎤ ⎡M x ′ ⎤ ⎢  ⎥ ⎢ ⎥ u t v t t M ( ) Ω ( ) Ω Δ ω ( r , ) Ω = − − − ( x′ y′ z′ ) M y ′ , 0 ⎢ y′ ⎥ ⎢ ⎥ ⎢  ⎥ ⎢ ⎣ M z ′ ⎥⎦ ⎢⎣ M z ′ ⎥⎦

(4)

⎡0 0 0 ⎤ ⎡ 0 0 1⎤ ⎡ 0 1 0⎤ Ω x ′ = ⎢0 0 −1⎥ , Ω y ′ = ⎢ 0 0 0 ⎥ , Ωz ′ = ⎢ −1 0 0 ⎥ , ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢⎣0 1 0 ⎥⎦ ⎢⎣ −1 0 0 ⎥⎦ ⎢⎣ 0 0 0 ⎥⎦

(7)

and we now assume that Δω0 = Δω (r , t ) = Δω (z , t ) as only gradient in z-direction is applied and it dominates the Δωoff and δω0 which are negligible. Eq. (6) clearly indicates that the excitation and inhomogeneities cause the magnetisation vector to rotate about an axis.1 If gradient fields are superimposed on the main static field then only the G r (t ) ⋅ r part of Δω0(r,t) is taken into account. Given that the gradient field is applied in the z-direction, the objective is to find controls, u(t) and v(t), that can drive the bulk magnetisation to the desired slice profile such that 2

J = ∫ M f (T f , u[0,T f ], v[0,T f ], z) − Md (z) dz

Here, ω1(t) is referred to as the Rabi frequency [32–34].

z

(8)

is minimised. In the above equation, Mf is the final state after the pulse has been turned off at t = Tf and Md represents the desired state.

Control of a spin system by the Fourier synthesis method For a short period of excitation (less than one time constant ( t  T 2 (r) ), it is possible to ignore relaxation terms [35] and approximate the magnetic resonance phenomenon in the classical rotating frame of reference with the following equation:

⎡M x ′ ⎤ ⎡ Δω0 (r , t ) v (t ) ⎤ ⎡M x ′ ⎤ 0 ⎢  ⎥ ⎢ M Δ ω ( , ) t u (t ) ⎥ ⎢M y ′ ⎥ , r = − 0 ′ y 0 ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢  ⎥ ⎢ −v (t ) 0 ⎥⎦ ⎢⎣ M z ′ ⎥⎦ −u (t ) ⎣⎢ M z ′ ⎦⎥ ⎣

(6)

in which

where

Δω0 (r , t ) = (ω0 − ωrf ) + γδ B 0 + γ G r (t ) ⋅ r = Δωoff + δω0 + Δω (r , t ).

Table 1 Description of the Bloch equation parameters.

Δω(r,t)

The behaviour of an ensemble of spins at a classical level in the presence of external magnetic fields may be described by the Bloch equation [31]. The Bloch equation in the classical rotating frame of reference whose transverse plane is rotating clock-wise at the Larmor frequency of the static field, ω0, is written as

95

Review of the Fourier synthesis method for slice selection in Cartesian coordinates The Fourier synthesis method may be used to solve the optimisation problem indicated by Eq. (8). In this technique during

(5) 1 exp(αΩ ) exp(βΩ ) x′ , y ′ , and exp(γΩ z ′ ) generate rotation matrices about x′, y′, and z′ axes, respectively.

96

B. Tahayori et al./Physica Medica 32 (2016) 94–103

series coefficients of that function. The number of the terms in Fourier series expansion, N, is determined based on the desired accuracy. The period of f(z) is considered to be 2π. Improving the time efficiency of the Fourier synthesis method using the Bloch equation in spherical coordinates We have demonstrated in our previous work [36,37] that the Bloch equation in spherical coordinates, where

M x ′ = M R ′ sin M φ ′ cos M θ ′,

(14a)

M y ′ = M R ′ sin M φ ′ sin M θ ′,

(14b)

M z ′ = M R ′ cos M φ ′

(14c)

Figure 1. The sequence of pulses for slice selective excitation.

the excitation period a mixture of gradient field and the excitation field is applied. This sequence is repeated several times (Fig. 1) to excite a desired slice profile.2 The excitation field in the laboratory frame of reference is applied along the x-axis, and therefore v(t) = 0. The sequence of pulses depicted in Fig. 1 causes a total rotation of bulk magnetisation as a result of each rectangular pulse excitation as well as the gradient fields. Thus, the net rotation will be k =N

U = ∏U k ,

for short pulses when MR′ = 1, can be written as

M θ ′ = −Δω + u (t )cot M φ ′ cos M θ ′ + v (t )cot M φ ′ sin M θ ′,

(15a)

M φ ′ = u (t )sin M θ ′ + v (t )cos M θ ′,

(15b)

with

M θ ′ (t = 0) = M θ0′ ,

(16a)

M φ ′ (t = 0) = M φ0′ .

(16b)

(9)

k =0

Here Mθ′ and Mφ′ are azimuth and colatitude components of the magnetisation vector, respectively. When u (t ) = v (t ) = 0 ,

where each Uk causes an incremental rotation

t

⎛ 1 ⎞ U k = exp(−α k Ω z′ )exp ⎜ − βk Ω x′ ⎟ exp(2α k Ω z′ ) ⎝ 2 ⎠ ⎛ 1 ⎞ exp ⎜ − βk Ω x′ ⎟ exp(−α k Ω z′ ), ⎝ 2 ⎠

(10)

Mφ ′ (t) = Mφ0′ .

(17b)

Moreover, when v(t) = 0 and Δω = 0, the solution to the dynamics described by Eq. (15) may be written as

τu

G(t) and u(t) at each step are τG and τu, respectively. In [26–28], it is shown that the net effect of the sequence of pulses represented in Fig. 1 can be approximated by k =N

U = ∏ U k ≈ exp(−f (z )Ω x ′ ),

(17a)

0

where αk = zAk, Ak = ∫ γ G (τ )dτ , and βk = ∫ u (τ )dτ . The durations of τG

Mθ ′ (t) = − ∫ Δω dτ + Mθ0′ ,

(11)

k =0

tan Mθ ′ = cos β (t) tan Mθ0′ + sin β (t)

cot Mφ0′ , cos Mθ0′

cos M φ ′ = cos β (t ) cos M φ0′ − sin β (t ) sin M φ0′

(18a)

sin M θ0′ ,

(18b)

t

where β (t ) = ∫ u (τ )dτ . The validity of the above equation can be in0

vestigated through substituting the above solutions to Eq. (15).

where k =N

k =N

k =0

k =0

f (z ) = ∑ βk cos α k = ∑ βk cos(Ak z ).

(12)

If the area under the gradient field, Ak, is selected such that Ak = k, then k =N

f (z ) = ∑ βk cos(kz ).

(13)

k =0

Thus, the proposed sequence of pulses causes the bulk magnetisation to rotate about the x′-axis with an angle f(z) which clearly is position dependent. Thus, to obtain a desired excited profile, fd(z), this function should be approximated by f(z). Eq. (13) states that βk coefficients are the Fourier Series coefficients of an even function. Thus for a given excited profile, fd(z), first we should define an even function, fde (z) , based on fd(z) and then calculate the Fourier

2 The Fourier synthesis method may be used to selectively excite a predetermined volume in the object. This result can be found in [26,27].

A novel sequence of pulses for slice selection Consider the novel sequence of pulses shown in Fig. 2. To calculate the azimuth and colatitude components at the end of this sequence, we adopt Eqs. (17) and (18). If the initial condition at the start of the process is

Ψ k− =2m ≡ Ψ k =2m (t = Tk 0 ) = [M θ−k′

M φ−k′ ] , T

(19)

then

Ψk (t = Tk1) = [ Mθk′ 1

Mφk′ 1 ] , T

(20)

where

tan M θk′ 1 = cos βk

tan(M θ−k′ ) + sin βk

cot M φ−k′ , cos M θ−k′

cos M φk′ 1 = cos βk

cos M φ−k′ − sin βk

sin M φ−k′

(21a)

sin M θ−k′ ,

(21b)

B. Tahayori et al./Physica Medica 32 (2016) 94–103

97

which results in4

⎛ a cos Mθ−k′ ⎞ Mθk′ 1  Mθ−k′ + tan−1 ⎜ k , ⎝ 1 + ak sin Mθ−k′ ⎟⎠

(27a)

M φk′ 1  M φ−k′ + βk

(27b)

sin M θ−k′ ,

where ak = βk cot M φ−k′ . It should be noted that the designed sequence always starts from U0 and therefore in the initial step, where k = 0 and as a result the gradient field is zero, the excitation field tips the bulk magnetisation immediately and during the sequence of pulses we will not be at M φ ′ = 0 . If the area under the gradient field is selected such that A1k = 2k + 1 and A2k = 2k + 3, we may write Figure 2. One step of the novel sequence to improve the time efficiency of the Fourier synthesis method.

Ψ k (t = Tk 2) = [M θk′ 2

M φk′ 2 ]

T T

⎡ ⎛ a cos M θ−k′ ⎞ ⎤ M φ−k′ + βk sin M θ−k′ ⎥ .  ⎢M θ−k′ − (2k + 1)z + tan−1 ⎜ k − ⎟ ⎣ ⎝ 1 + ak sin M θk′ ⎠ ⎦ Ψ k (t = Tk 2) = [M θk′ 2

M φk′ 2 ]

T T

⎛ a cos M θ−k′ ⎞ ⎡ ⎤ M φ−k′ + βk sin M θ−k′ ⎥ (28)  ⎢M θ−k′ − (2k + 1)z + tan−1 ⎜ k − ⎟ ⎝ 1 + ak sin M θk′ ⎠ ⎣ ⎦

in which βk = ∫ u(τ )dτ = ukτ uk . Moreover, it can be shown that τ uk

Ψk (t = Tk 2) = [ Mθk′ 2

Mφk′ 2 ] = [ Mθk′ 1 − A1k z Mφk′ 1 ] , T

T

(22)

where the magnitude of the gradient field is chosen to be

A 1k = ∫

τG k +1

γ G (τ )dτ .

At t = Tk 3 the components of the magnetisation will be

Ψk (t = Tk 3) = [ Mθk′ 3

Mφk′ 3 ] , T

(23)

where

tan Mθk′ 3 = cos βk+1

cot Mφk′ 2 tan(Mθk′ 2) + sin βk+1 , cos Mθk′ 2

cos M φk′ 3 = cos βk +1 cos M φk′ 2 − sin βk +1 sin M φk′ 2 In the above equation βk+1 = ∫

τ uk +1

sin M θk′ 2 .

Mφk′ 4 ] = [ Mθk′ 3 + A2k z Mφk′ 3 ] , T

(25)

in which the magnitude of the gradient field is chosen such that

A 2k = ∫

τGk + 2

γ G (τ )dτ .

⎞ −1 ⎛ ak +1 cos M θk′ 2 ⎟ + tan ⎜ ⎝ 1 + ak +1 sin M θk′ 2 ⎠

⎞ ⎟, ⎠

⎛ a cos Mθk′ 2 ⎞ Mθk′ 3  Mθk′ 2 + tan−1 ⎜ k+1 ⎝ 1 + ak+1 sin Mθk′ 2 ⎟⎠ ⎛ a cos Mθ−k′ ⎞  Mθ−k′ − (2k + 1) z + tan−1 ⎜ k ⎝ 1 + ak sin Mθ−k′ ⎟⎠ ⎛ a cos Mθk′ 2 ⎞ + tan−1 ⎜ k+1 ⎝ 1 + ak+1 sin Mθk′ 2 ⎟⎠

(24b)

u(τ )dτ = uk+1τ uk+1 . The final value

T

⎛ a cos M θk′ 2 ⎞ M θk′ 3  M θk′ 2 + tan−1 ⎜ k +1 ⎟ ⎝ 1 + ak +1 sin M θk′ 2 ⎠ ⎛ a cos M θ−k′  M θ−k′ − (2k + 1)z + tan−1 ⎜ k − ⎝ 1 + ak sin M θk′

(24a)

of the azimuth and colatitude components at step k is calculated to be

Ψk+ = Ψk (t = Tk 4 ) = [ Mθk′ 4

Similarly, Eq. (24) leads to

M φk′ 3  M φk′ 2 + βk +1 sin M θk′ 2  M φ−k′ + βk

(29a)

sin M θ−k′ + βk +1 sin M θk′ 2 , (29b)

where ak +1 = βk +1 cot M φk 2 . The components of the magnetisation at the end of the sequence of pulses represented by Fig. 2 will be

⎡M θ ′ ⎤ Ψ k+ = ⎢ k 4 ⎥ ⎣M φk′ 4 ⎦ ⎛ a cos M θ−k′ ⎞ ⎡ − ⎛ a cos M θk′ 2 M θk′ + 2z + tan−1 ⎜ k + tan−1 ⎜ k +1 ⎢ − ⎟ ⎢ ⎝ 1 + ak sin M θk′ ⎠ ⎝ 1 + ak +1 sin M θk′ 2 ⎢⎣ M φ−k′ + βk sin M θ−k′ + βk +1 sin M θk′ 2

⎞⎤ ⎟ ⎥ (30) ⎠⎥ . ⎥⎦

Approximations Calculating the final values of magnetisation components For small values of βk, Eq. (21) may be approximated by3

tan Mθk′ 1 = tan(Mθ−k′ ) + βk cos M φk′ 1 = cos M φ−k′ − βk

cot Mφ−k′ , cos Mθ−k′ sin M φ−k′

(26a)

sin M θ−k′ ,

Consider Fig. 3 which consists of the application of m = N/2 consecutive sequences of pulses represented by Fig. 2 to a spin system.5 In the Appendix, we have proven that if the initial condition of the system is

(26b) a a cos x , then it can be shown that y = mπ + x + tan−1 cos x 1 + a sin x where m is an integer. Moreover, if cos(y ) = cos(x ) , then y = 2mπ ± x. Without loss of generality we have assumed that m = 0. 5 N is an even non-negative integer. 4

3 All approximations made in this section can be replaced by an equality if O(β 2 ) k is added to the right hand side which justifies the mathematical induction used to prove the validity of Fourier synthesis method.

If tan y = tan x +

98

B. Tahayori et al./Physica Medica 32 (2016) 94–103

Figure 3. Application of m = N/2 consecutive steps of the sequence of pulses represented in Fig. 2.

Ψ 0− ≡ [M θ−0′

M φ−0′ ] = [π 2 M φ−0′ ] , T

T

(31)

the objective is to show that at the end of the sequence of pulses the azimuth and colatitude components will be

⎡ ⎢ Ψ ⎢ ⎢ Mφ−0′ ⎣ ⎡ ⎢ =⎢ ⎢⎣ Mφ−0′ + N

π + (N + 2) z 2

⎤ ⎥ ⎥ k=N + ∑ k=0,2,… βk cos kz + βk+1 cos(k + 1) z ⎦⎥ π ⎤ + (N + 2) z ⎥ 2 ⎥. k=N + ∑ βk cos kz ⎥⎦

(32)

k =0

Novel sequence of pulses for Fourier synthesis method From Eq. (32), the colatitude component after application of the sequence of pulses shown in Fig. 3 will be k=N

Mφ′ (z) = Mφ−0′ + ∑ βk cos kz ,

(33)

k =0

which can be interpreted that M φ′ (z ) can be set equal to a truncated Fourier series of N terms where each term is small (i.e. βk  1 ). As a result, to obtain a desired excitation profile, M φd′ (z ) , we have to define an even function Mφde′ (z) and calculate its Fourier expansion coefficients which provides us the magnitude of the excitation pulses in Fig. 4. To provide a uniform excitation, the azimuth component must be independent of position. As can be seen from Eq. (32), the azimuth component is a function of z as well. However, it is possible to remove this space dependency by applying a gradient field with a proper magnitude such that

∫τ

g

γ G(τ )dτ = N + 2 . Therefore, the pulse

sequence that can selectively excite a slice consists of m = N/2 sequences of pulses shown in Fig. 2 followed by a gradient field such that it removes the space dependency of the azimuth component. It is obvious from Fig. 4 that the last two steps can be combined and a gradient whose area is N + 1 may be applied in the negative z′ direction. Remarks on the proposed sequence of pulses If the same number of terms is used for the original sequence depicted in Fig. 1 and the novel proposed sequence of pulses shown in Fig. 4, it is expected that they have similar performance. In the Results section, we will investigate this qualitatively as well as quantitatively. The interesting point is that for each coefficient, βk, of the original sequence of pulses, four gradient pulses each with area k are required. However, in the proposed sequence, for every βk, only one gradient pulse with area 2k + 1 is required. Therefore, the time required for rotations made by gradient field for the novel sequence of pulses is roughly half of that for the original sequence. To be more accurate, the overall required rotation in radians both in positive and negative directions for the original sequence of pulses, if N terms in Fourier series expansion are kept, will be

RGOriginal = 4(1 + 2 + 3 + …+ N) z = 2N (N + 1) z .

(34)

However, this value for the novel sequence of pulses is

RGNovel = (1 + 3 + 5 + … + (2N − 1) + (2N + 1) + (N + 1)) z = (N + 1)(N + 2) z 

1 RGOriginal . 2

(35)

Since the time required to make these rotations is linearly proportional to the rotation value in radian, then the total time will be halved as well. It is important to note that the time required for the excitation field to induce the designed rotation is the same for both sequences if the same number of terms in Fourier series is kept. The overall rotation required by the excitation field is the summation of all Fourier coefficients which is much smaller than the one required by gradient fields. Therefore, the dominant factor in calculating the duration of N steps of the sequence of pulses is the required gradient field time. Thus, the proposed sequence of pulses has an improved time efficiency compared to the original one. It should be noted that a similar sequence of pulses can be used to suppress the RF field inhomogeneity. The performance of the novel sequence of pulses for RF inhomogeneity suppression is similar to that of the original sequence, nonetheless, the novel sequence will be more time efficient. Results

Figure 4. The novel sequence of pulses for slice selective excitation for Fourier synthesis method.

To obtain a better insight into the Fourier synthesis method, we have simulated the original sequence of pulses as well as the novel one proposed in this paper. It should be noted that all simulations presented in this section are performed by using the original

B. Tahayori et al./Physica Medica 32 (2016) 94–103

99

Figure 5. Simulation results to excite a slab of 2 mm using Fourier synthesis method: (a) original sequence of pulses with N = 50, (b) original sequence of pulses with N = 100, (c) novel sequence of pulses with N = 50 and (d) novel sequence of pulses with N = 100. The desired rotation about the x′-axis is set to π/2.

equations and not the approximate results which validate our approximations. Applications of the Fourier synthesis method for slice selection In this section a few applications of the Fourier synthesis method will be presented to show how the original Fourier synthesis sequence of pulses controls the spin system for slice selection. Simulation results will be presented both in Cartesian and spherical coordinates. Application 1: It is desired to excite a slice of an object in the following form:

⎧α d fd (z) = ⎨ ⎩0

z1 < z < z2 otherwise

(36)

where αd determines the rotation angle about the x′-axis and z2 > z1 > 0. We define fde (z) to be

⎧α d fde (z) = ⎨ ⎩0

z1 < z < z2 . otherwise

(37)

By Fourier series expansion the values of βk are calculated as follows:

βk =

2α d sin kz0, k = 1, 2, 3, …. kπ

(40b)

Figure 5a and b represents the simulation result of Application 1 for a slab of thickness 2 mm in Cartesian coordinates using the original sequence of pulses, when 50 and 100 Fourier series coefficients are used, respectively. The desired rotation angle is set to be αd = π/2. The simulation result for the novel sequence of pulses is shown in Fig. 5c and d for 50 and 100 Fourier series coefficients, respectively. Application 2: It is desired to excite two slices of an object in the following form:

⎧αd ⎪ f d (z ) = ⎨αd ⎪0 ⎩

z1 < z < z 2 z3 < z < z4 otherwise

(41)

where α d is the desired rotation angle about the x′-axis and z 4 > z 3 > z 2 > z 1 > 0 . In this case we have

⎧α d ⎪ fde (z) = ⎨α d ⎪0 ⎩

z1 < z < z2 z3 < z < z 4 otherwise

(42)

β0 =

αd (z 2 − z 1), π

(38a)

By Fourier series expansion the values of βk are calculated as follows:

βk =

2αd (sin kz 2 − sin kz 1), k = 1, 2, 3, …. kπ

(38b)

β0 =

αd (z 2 − z 1 + z 4 − z 3), π

(43a)

If the desired slice has symmetry about z = 0, its profile may be written as

βk =

2αd (sin kz 2 − sin kz 1 + sin kz 4 − sin kz 3) k = 1, 2, 3, …. kπ

(43b)

⎧α d fd (z) = ⎨ ⎩0

Figure 6 represents the simulation result of Application 2 in Cartesian coordinates when z1 = 2 mm, z2 = 4 mm, z3 = 8 mm, and z4 = 10 mm. In this simulation 100 terms are taken into account in Fourier series expansion. The desired rotation angle is set to be αd = π/ 2. This simulation clearly shows that it is possible to excite more than one slice of the object in the same amount of time required to excite one slice. Moreover, this goal can be achieved by applying a single frequency RF pulse.

z < z0 otherwise,

(39)

then f de ≡ f d and, therefore, z2 = z0 and z1 = 0. Thus, the Fourier series coefficients are

β0 =

αd z 0, π

(40a)

100

B. Tahayori et al./Physica Medica 32 (2016) 94–103

Table 2 Performance of the original Fourier synthesis sequence of pulses versus the improved Fourier synthesis sequence of pulses. The root mean squared error of the magnetisation vector, M, is presented. Original sequence of pulses N

10 20 50 100

Inside the slice

Outside the slice

Overall

Total

RMSE

RMSE

RMSE

Time (ms)

0.281 0.139 0.103 0.087

0.311 0.206 0.132 0.091

0.419 0.248 0.168 0.126

1.4 4.4 24.6 96.3

Inside the slice

Outside the slice

Overall

Total

RMSE

RMSE

RMSE

Time (ms)

0.263 0.132 0.099 0.085

0.321 0.209 0.136 0.095

0.415 0.247 0.169 0.128

1 2.6 13.0 49.3

Figure 6. Simulation results of the novel Fourier synthesis sequence of pulses in Cartesian coordinates to simultaneously excite two slabs of the object, z1 = 2 mm, z2 = 4 mm, z3 = 8 mm, and z4 = 10 mm with N = 100 and αd = π/2.

Novel sequence of pulses

Application 3: It is desired to excite a trapezoidal slice with the following form:

10 20 50 100

⎧α d ⎪ (z − z) ⎪ αd fd (z) = ⎨ 1 ⎪ z1 − z0 ⎩⎪0

0 < z < z0 z0 < z < z1,

(44)

otherwise

the Fourier series are calculated to be

β0 = βk =

N

αd z −z z 0 + 1 0 αd , π 2π

(45a)

2α d 2α (cos kz1 − cos kz0 + kz1 sin kz0 − kz0 sin kz0 ) sin kz0 + d kπ k 2π (z1 − z0 ) k = 1, 2, …. (45b)

Since a trapezoidal function is smoother than a rectangular function the convergence rate is quicker. The k2 term in the denominator of Eq. (6) indicates this faster convergence rate. The simulation result for 30 terms in Fourier series expansion for a slice with z0 = 1 mm and z0 = 2 mm is shown in Fig. 7. As can be observed the performance of the method for transverse magnetisation Mx ′y ′ is as good as the rectangular excitation with 50 terms.

Comparison of the original and novel sequences for Fourier synthesis method To make a comparison between the performance of the original and the novel sequences of pulses for slice selection, we compare their results with an ideal slice which has a rectangular excited profile. Moreover, we compare the time efficiency of these sequences versus each other quantitatively.

Comparing the performance of both sequences of pulses To interpret the performance of Fourier synthesis method slice selectivity for symmetric case when z0 = 1 mm, we compare their result with an ideal rectangular profile. Figure 8 shows the squared error between the ideal and the excited slice profiles shown in Fig. 5b and d. As we expected the performance of the Fourier synthesis method reduces at the edges of the desired slice.6 The qualitative comparison of Fig. 8a and b states that the performance of both sequences of pulses is almost the same. To compare the results quantitatively we compare them in terms of the root mean squared error (RMSE). Table 2 represents the root mean squared error of both methods for different number of terms in Fourier series expansion. Comparison of the values in this table states that the performance of the novel sequence of pulses in sense of RMSE is as good as the original sequence of pulses. Comparing the duration of both sequences of pulses Consider a gradient field of strength 1000 mT/m which is typical gradient field strength for an animal scanner. Assuming the rise time and the fall time are 1 μs each, this gradient field rotates the bulk magnetisation for 1 rad/mm in approximately 4.75 μs for the resonance frequency of hydrogen atom.7 A 1 mT excitation field causes a 1 rad rotation in approximately 3.7 μs. If we assume that the duration of each uk in the sequence of pulses is set to 4 μs and 50 terms are kept in the Fourier series expansion, the overall time required by the original sequence of pulses is calculated to be 24.6 ms. However, for the proposed sequence of pulses the duration of the sequence will be 13 ms. As expected, the time efficiency of the method is approximately double while maintaining the same performance. Table 2 shows the total time required by both sequences of pulses for different number of terms kept in the Fourier series expansion. RF field inhomogeneities and slice selectivity RF field inhomogeneities can affect the performance of any slice selective excitation, specifically at high field magnets. In the presence of the RF field inhomogeneities, the excitation pattern in the Bloch equation will be λu(t) where λ ∈[1 − ε1, 1 + ε 2] , 0 < ε1, ε 2 < 1.

Figure 7. Simulation results of the novel Fourier synthesis sequence of pulses in Cartesian coordinates to excite a trapezoidal slab of 2 mm, z0 = 1 mm and z1 = 2 mm. In Fourier expansion of the desired profile we have considered N = 30 terms. The desired rotation about the x′-axis is chosen to be π/2.

6 7

All slice selective pulse design algorithms experience the same problem. γ for hydrogen is 42.58 MHz/T. Given that 2πγ × 1000 × 10−3 × 1× 10−3 (τ G − 1) = 1,

τG is found to be 37.5 μs.

B. Tahayori et al./Physica Medica 32 (2016) 94–103

(a)

101

(b)

Figure 8. Squared error of the Fourier synthesis method result depicted in Fig. 5 in comparison with an ideal slice profile when N = 100: (a) original sequence of pulses and (b) novel sequence of pulses.

To investigate the effect of RF field inhomogeneities on the performance of the original and novel sequences of pulses, we have simulated the Fourier synthesis method for both sequences of pulses with N = 100 when λ = 0.8 and λ = 1.2. The results are demonstrated by Fig. 9. Field inhomogeneities, assumed to be constant in one slice, change the desired angle of rotation linearly. This can be inferred from the Fourier series expansion of a pulse that directly depends on the angle of rotation, αd, e.g. see Eq. (38). The overall RMSE (logarithmic scale) as a function of λ for both sequences of pulses is presented in Fig. 10 when N = 100. The results illustrate that the performance of the slice selection algorithm decreases by the RF field inhomogeneities. However, both sequences show similar behaviour as the RF field inhomogeneity increases. The overall RMSE is calculated for the magnetisation vector, M. A figure of merit for gradient complexity To measure the gradient complexity, we have introduced the following figure of merit:

Tp Gradient Complexity 

1 Tp Tp dt (∂G ∂t)2 Tp ∫ dt (∂G ∂t)2 (46) Tp ∫0 0 = , Gmax Gmax

where G max is the maximum magnitude of the applied gradient and Tp is the overall time for a sequence of pulses applied to select a

slice. The time derivative of the gradient field with regard to time is an indicator of how fast the gradient should switch. Figure 11 shows the natural logarithm of normalised gradient complexity versus RMSE for slice selection. This figure clearly shows that to selectively excite more uniform slices a more complicated gradient field is required for both the original and the improved Fourier synthesis methods. However, for a desired uniformity in the selected slice, the improved Fourier synthesis method requires a less complicated Gradient field.

Discussion The sequence of pulses presented in this paper to selectively excite a thin slice of the object requires fast switching gradient fields. Specifically, as the number of terms used in the Fourier series expansion of the desired slice profile increases, the frequency of the switching increases as well. This fast switching is required as the derivation was based on the assumption that the overall duration of the slice selection is sufficiently short to neglect the effect of relaxation time constants. Therefore, there is a trade-off between slice selection accuracy and the switching frequency of the gradient fields. It should be noted that with current technology there are limits to the gradient field switching frequencies which may restrict the application of the proposed method. If gradient fields can be switched more quickly in the future, one will be able to include sufficient terms

(a)

(b)

(c)

(d)

Figure 9. Simulation results to excite a slab of 2 mm using Fourier synthesis method in the presence of RF field inhomogeneities with N = 100 for (a) original sequence of pulses with λ = 0.8, (b) original sequence of pulses with λ = 1.2, (c) novel sequence of pulses with λ = 0.8 and (d) novel sequence of pulses with λ = 1.2. The desired rotation about the x′-axis is set to π/2.

102

B. Tahayori et al./Physica Medica 32 (2016) 94–103

Figure 10. The RMSE for both sequences of pulses in the presence of RF field inhomogeneities characterised by λ. N = 100 for both sequences of pulses.

Logarithm of Normalised Gradient Complexity

1.1 Oroginal F.S. Improved F.S

1 0.9 0.8 0.7

magnetic resonance. Fourier synthesis method that has been introduced recently uses a combination of the RF and gradient fields to selectively excite the spins in an object. In this paper, we reviewed the Fourier synthesis method for slice selection in MRI. In this method for slice selection, the gradient fields as well as the RF field are switched to selectively excite a thin layer of spins in the object under imaging. The method can also be used to suppress the RF field inhomogeneities. We proved the validity of the method in Cartesian coordinates as well as the spherical coordinates. The Bloch equations in Cartesian coordinates and spherical coordinates were used to provide a proof for this method of slice selection. Based on the spherical coordinates proof, we proposed a novel sequence of pulses that uses the Fourier synthesis method to selectively excite a slice of the object under imaging. The major advantage of this novel sequence of pulses is that it excites a slice in half of the time required for the original sequence of pulses. We illustrated the proof through numerical simulation and demonstrated how the Fourier synthesis method can be used for slice selection. We introduced a gradient complexity measure to compare the performance of the original and the proposed sequences of pulses for slice selection. Simulation results demonstrate that the original Fourier synthesis method has a higher gradient complexity for the same level of uniformity in the excited slice.

0.6 0.5 0.1

Acknowledgment

0.15

0.2

0.25

0.3

0.35

0.4

0.45

RMSE for slice selection Figure 11. Natural logarithm of normalised gradient complexity versus RMSE for slice selection.

in the Fourier series expansion which will result in improving the practical performance of the proposed sequence of pulses. It is worth mentioning that in this approach the gradient fields do not need to have a rectangular shape. The main factor is the area under the applied gradient field. Therefore, smoother gradients, e.g. sinusoidal gradients, may be applied which makes the results obtained through this method more achievable. Numerical simulations illustrate that the result of RF field inhomogeneities is a poor slice selectivity using the Fourier synthesis method. The desired angle of rotation alters linearly with RF field inhomogeneities and, therefore, the quality of the slice selection decreases. However, both sequences of pulses show similar behaviour as a function of the RF field inhomogeneities. The derivation of the Fourier synthesis method for both sequences of pulses is based on the assumption that the RF field as well as the gradient field rotate the bulk magnetisation while keeping its magnitude. This assumption is only true when the relaxation time constants in the Bloch equation can be ignored. However, for tissue types with very low T2, e.g. cortical bone or skeletal muscle, the pulse duration is comparable to T2 and thus the basic assumption of the Fourier synthesis method is no longer valid. Therefore, the slice selection performance is expected to be poor. A potential problem with both the original and the improved Fourier synthesis methods is the complexity of the associated gradient fields which can limit practical applications compared to the popular SLR technique. However, the slice selection uniformity of the original and the improved Fourier synthesis methods can be more selective than that of SLR when a fast switching gradient can be applied, which captures the potential benefit of this method. Conclusions The design of slice selective pulses for magnetic resonance imaging is still an active area of research in the field of nuclear

This work was supported by the National ICT Australia Victoria Research Laboratory. The first author received Postgraduate Overseas Research Experience Scholarships (PORES) from the University of Melbourne to spend four months at School of Engineering and Applied Mathematics of Harvard University. The authors would like to thank the reviewers for their insightful comments that considerably improved the quality of the paper. Appendix Proof of the Fourier synthesis method for the novel sequence of pulses in spherical coordinates using mathematical induction We use mathematical induction to prove the validity of the solution given by Eq. (32). Since N is an even number, we will show the validity of the statement for P0 and Pn when we assume that Pn−2 is true. P0: For the initial step we investigate the validity of

π Ψ 0+  ⎡⎢ + 2z ⎣2

T

M φ−0′ + β0 + β1 cos z ⎤⎥ . ⎦

(A1)

By the assumption made for the initial conditions of the system, M θ−0′ = π 2 , and using Eq. (28) when k = 0 leads to

M θ02 ′ 

π −z. 2

(A2)

By adopting the results expressed in Eq. (30), it can be shown that

π ⎡ ⎤ + 2z + O(β12) ⎢ ⎥ ⎡ Mθ04 ′ ⎤ 2 Ψ =⎢ ⎢ ⎥ ⎥ ′ ⎦ ⎣ Mφ04 ⎢ M − + β sin π + β sin ⎛ π − z ⎞ ⎥ φ0′ 0 1 ⎜ ⎟ ⎢⎣ ⎝ 2 ⎠ ⎥⎦ 2 π ⎡ ⎤ + 2z ⎥. ⎢ 2 ⎢ − ⎥ β β M cos z + + 0 1 ⎣ φ0′ ⎦ + 0

(A3)

B. Tahayori et al./Physica Medica 32 (2016) 94–103

The above results prove the validity of the statement for the initial step. Pn−2 : If the statement is true for step n − 2, the azimuth and colatitude components at the end of this step can be written as T

k =n− 2 ⎡π ⎤ Ψn+−2  ⎢ + nz Mφ−0′ + ∑ k=0,2,... βk cos kz + βk+1 cos(k + 1) z ⎥ . ⎣2 ⎦

(A4)

The final values at step n − 2 are the initial values of step n, therefore we may write T

k =n− 2 ⎡π ⎤ Ψn−  ⎢ + nz Mφ−0′ + ∑ k=0,2,… βk cos kz + βk+1 cos(k + 1) z ⎥ . ⎣2 ⎦

(A5)

Pn: By replacing the initial conditions indicated by Eq. (A5) in Eq. (28) when k = n, we obtain

M θn′ 2 

π π − (n + 1)z + O (βn2)  − (n + 1)z . 2 2

(A6)

Using Eq. (30) results are shown in which π ⎡ ⎤ + (n + 2)z + O (βn2) ⎢ ⎥ 2 Ψ n+  ⎢ ⎥ k =n −2 − ⎢⎣M φ0′ + ∑ k =0,2,... ( βk cos kz + βk +1 cos(k + 1)z ) + βn cos nz + βn +1 cos(n + 1)z ⎥⎦ π π ⎡ ⎤ ⎡ ⎤ + (n + 2)z + (n + 2)z ⎢ ⎥ ⎢ ⎥ 2 2 =⎢ ⎢ (A7) ⎥ ⎥. k =n k =N ⎢⎣M φ−0′ + ∑ k =0,2,... βk cos kz + βk +1 cos(k + 1)z ⎥⎦ ⎢⎣M φ−0′ + ∑ βk cos kz ⎥⎦ k =0

Thus, the statement is true for n and the proof is complete. Specifically, if M φ−0′ = 0 at the start of the pulse sequence, then the final position of the bulk magnetisation will be

π ⎡ ⎤ ⎡ π ⎤ + (N + 2) z ⎢ ⎥ ⎢ 2 + (N + 2) z ⎥ 2 Ψ +N  ⎢ = ⎥ ⎢ k=N ⎥. k=N ⎢∑ ⎥ ⎢ ⎥⎦ cos kz + cos( k + ) z 1 β β kz cos β k +1 ⎣ k=0,2,… k ⎦ ⎣ ∑ k =0 k

(A8)

References [1] Westbrook C, Roth C, Talbot J. MRI in practice. 4th ed. Blackwell; 2011. [2] Liang ZP, Lauterbur P. Principles of magnetic resonance imaging, a signal processing perspective. Piscataway, NJ, USA: IEEE Press; 2000. [3] Garroway A, Grannell P, Mansfield P. Image formation in NMR by a selective irradiative process. J Phys 1974;7:457–62. [4] Balchandani P, Pauly J, Spielman D. Slice-selective tunable-flip adiabatic low peak-power excitation (stable) pulse. Magn Reson Med 2008;59:1072–8. [5] Moore J, Jankiewicz M, Zeng H, Anderson AW, Gore JC. Composite rf pulses for -insensitive volume excitation at 7 Tesla. J Magn Reson 2010;205(1):50–62. [6] Finsterbusch J, Busch MG. Segmented 2D-selective RF excitations with weighted averaging and flip angle adaptation for MR spectroscopy of irregularly shaped voxel. Magn Reson Med 2011;66(2):333–40. [7] Moore J, Jankiewicz M, Anderson AW, Gore JC. Slice-selective excitation with -insensitive composite pulses. J Magn Reson 2012;214:200–11. [8] Park B, Sung K, McGarrity J, Oh SK, Wang Z, Collins CM. Slice-selective transmit array pulses for improvement in excitation uniformity and reduction of sar. J Electromagn Anal Appl 2013;5:205–12.

103

[9] Balchandani P, Glover G, Pauly J, Spielman D. Improved slice-selective adiabatic excitation. Magn Reson Med 2014;71:75–82. [10] Hoult D. The solution of the Bloch equations in the presence of a varying B1 field, an approach to selective pulse analysis. J Magn Reson 1979;35:69–86. [11] Mansfield P, Maudsley A, Morris P, Pykett I. Selective pulses in NMR imaging, a reply to criticism. J Magn Reson 1979;33:261–74. [12] Conolly S, Macovski A. Selective excitation via optimal control theory. Proceedings of the 4th Meeting of the Society of Magnetic Resonance in Medicine, p. 958, 1985. [13] Locher P. Computer simulation of selective excitation in NMR imaging. Philos Trans R Soc Lond B Biol Sci 1980;289:537–42. [14] Loeffler W, Opplet A, Faul D. Computer simulations of slice selection in NMR imaging. Proceedings of the 2nd Meeting of the Society of Magnetic Resonance in Medicine, p. 196–7, 1983. [15] Hinshaw WS, Lent AH. An introduction to NMR imaging, from the Bloch equation to the imaging equation. Proc IEEE 1983;71:338–50. [16] Conolly S, Nishimura D, Macovski A. Optimal control solutions to the magnetic resonance selective excitation. IEEE Trans Med Imaging 1986;MI-5:106–15. [17] Pauly J, Le Roux P, Nishimura D, Macovski A. Parameter relations for the Shinnar-Le Roux selective excitation pulse design algorithm. IEEE Trans Med Imaging 1991;10:33–65. [18] Ulloa JL, Guarini M, Guesalaga A, Irarrazaval P. Chebyshev series for designing RF pulses employing an optimal control approach. IEEE Trans Med Imaging 2004;23:1445–52. [19] Xu D, King KF, Zhu Y, McKinnon GC, Liang ZP. Designing multichannel, multidimensional, arbitrary flip angle RF pulses using an optimal control approach. Magn Reson Med 2008;59:547–60. [20] Li JS. Control of inhomogeneous ensembles [Ph.D. thesis]. Harvard University; 2006. [21] Li J, Khaneja N. Control of inhomogeneous quantum ensembles. Phys Rev A 2006;73:030302. [22] Ikinomidou V, Sergiadis G. Improved Shinnar-Le Roux algorithm. J Magn Reson 2000;143:30–4. [23] Raddi A, Klose U. Optimized Shinnar-Le Roux RF 180 pulses in fast spin-echo measurements. J Magn Reson Imaging 1999;9:613–20. [24] Lee KJ. General parameter relations for the Shinnar-Le Roux pulse design algorithm. J Magn Reson 2007;186:252–8. [25] Balchandani P, Glover G, Pauly J, Spielman D. Improved slice-selective adiabatic excitation. Magn Reson Med 2014;71(1):75–82. doi:10.1002/mrm.24630. [26] Pryor B, Khaneja N. Fourier decompositions and pulse sequence design algorithms for nuclear magnetic resonance in inhomogeneous fields. J Chem Phys 2006;125:19411. [27] Pryor B. Fourier synthesis methods for identification and control of ensembles [Ph.D. thesis]. Harvard University; 2007. [28] Pryor B, Khaneja N. Fourier methods for control of inhomogeneous quantum systems. In Proceedings of the 46th IEEE Conference on Decision and Control, p. 6340–5, Dec. 2007. [29] Arthanari H, Frueh D, Wagner G, Pryor B, Khaneja N. Fourier synthesis techniques for NMR spectroscopy in inhomogeneous fields. J Chem Phys 2008;128:214503. [30] Owrutsky P, Khaneja N. Control of inhomogeneous ensembles on the Bloch sphere. Phys Rev A 2012;86:022315. [31] Bloch F. Nuclear induction. Phys Rev 1946;70:460–74. [32] Allen L, Eberly JH. Optical resonance and two-level atoms. 3rd ed. New York, USA: Wiley Interscience; 1975. [33] Layton KJ, Tahayori B, Mareels IMY, Farrell PM, Johnston LA. Rabi resonance in spin systems: theory and experiment. J Magn Reson 2014;242:136–42. [34] Tahayori B, Johnston LA, Layton KJ, Farrell PM, Mareels IMY. Solving the Bloch equation with periodic excitation using harmonic balancing: application to Rabi modulated excitation. IEEE Trans Med Imaging 2015;34(10):2118–30. [35] Haacke EM, Brown RW, Thompson MR, Venkatsen R. Magnetic resonance imaging, physical principles and sequence design. Wiley-Liss; 1999. [36] Tahayori B, Johnston LA, Mareels IMY, Farrell PM. Novel insight into magnetic resonance through a spherical coordinate framework for the Bloch equation. In Proceedings of the SPIE Conference on Medical Imaging, vol. 7258, p. 72580Y1–12, 2009. [37] Tahayori B. Magnetic resonance imaging and the Bloch equation: a reappraisal [Ph.D. thesis]. The University of Melbourne; 2010.