Basics of Chemical Exchange Saturation Transfer (CEST) Magnetic Resonance Imaging

2.1. Two-pool chemical exchange model

A two-pool chemical exchange model is illustrated in Figure 1 . A and B in Figure 1 represent the pools of bulk water protons and labile solute protons, respectively. The time-dependent Bloch-McConnell equations for the two-pool chemical exchange model in CEST or APT MRI are expressed as [17, 18].

where superscripts a and b show the parameters in pool A and pool B, respectively. For example, Mxat, Myat, and Mzatare the x, y, and z magnetization components in pool A at time t, respectively. R1aand R2aare the reciprocals of the longitudinal (T1a) and transverse relaxation times (T2a), that is, the longitudinal and transverse relaxation rates in pool A, respectively. kab and kba denote the exchange rate from spins in pool A to those in pool B and that from spins in pool B to those in pool A, respectively ( Figure 1 ). M0aand M0bare the thermal equilibrium z magnetization components in pool A and pool B, respectively. Δωa  = ωa  − ω and Δωb  = ωb  − ω, where ωa , ωb , and ω denote the Larmor frequencies in pool A and pool B, and the frequency of the RF-pulse irradiation, respectively. ω1xand ω1yare the x and y components of the amplitude of the RF-pulse irradiation (ω ₁), respectively. Note that ω ₁ = γB ₁, where γ and B ₁ are the gyromagnetic ratio (γ/2π = 42.58 MH_z/T) and RF power, respectively. When the RF pulse is applied along an angle φ from the x-axis of the rotating frame as illustrated in Figure 2 , ω1xand ω1yare represented by ω1x = ω ₁ cos ϕ and ω1y = ω ₁ sin ϕ, respectively. When the RF pulse is applied along the x-axis of the rotating frame, ω1xand ω1ybecome ω ₁ and 0, respectively.

The differential equations given by Eq. (1) can be combined into one vector equation (homogeneous linear differential equation) [18]:

where

and

T in Eq. (3) denotes the matrix transpose.

For simplicity, we assume that the RF pulse is applied along the x-axis of the rotating frame, that is, ϕ = 0. According to Koss et al. [19], the matrix A( ω ,  ω ₁, 0) can be given by.

where E is the evolution matrix and C is the constant-term matrix. Furthermore, E is given by.

In the case of A given by Eq. (4), R is reduced to.

where

and

K in Eq. (6) is given by

where I is a 3-by-3 identity matrix and ⊗ denotes the Kronecker tensor product. C in Eq. (5) is given by.

The solution of Eq. (2) with ϕ being 0 can be given by [18].

where t represents the so-called saturation time and M(0) is the matrix of initial values at t = 0. e ^{A( ω ,  ω ₁, 0)t} is the matrix exponential.

It should be noted that mass balance imposes the following relationship between the exchange rates (kab and kba ) of pool A and pool B [17]:

and

2.2. Three-pool chemical exchange model

Figure 3 illustrates a three-pool chemical exchange model in which pool a represents the bulk water pool. In this case, R and K are given by [19].

and

respectively. R ^c in Eq. (15) is given by Eq. (8) in which the subscript a and superscript a are replaced by c. C is given by.

The solutions of other multi-pool chemical exchange models such as an hour-pool chemical exchange model are described in Ref. [20].

2.3. Calculation of Z-spectrum, MTR_asym, and PTR

The CEST effect has usually been analyzed using the so-called Z-spectrum [18]. The Z-spectrum is given by the following equation:

where MzaΔωoffis the z magnetization component of bulk water protons (pool A) at Δωoff . Note that Δωoff  =  − Δωa .

The magnetization transfer asymmetry (MTR_asym) analysis has been performed using the following equation [18]:

Instead of MTR_asym, the following equation for proton transfer ratio (PTR) has also been used for analyzing the CEST effect [18]:

where Mza−Δωoffdenotes the z magnetization component of pool A at the opposite side of the water resonance (Δωoff ).

Figure 4(a) shows Z-spectra as a function of offset frequency (Δωoff ) for various saturation times (0.5, 1, 2, 5, and 10 s) in the two-pool chemical exchange model ( Figure 1 ). Figure 4(b) shows Z-spectra as a function of Δωoff for various ω ₁ values (25, 50, 100, 150, and 200 Hz). It should be noted that because B ₁ = ω ₁/γ, ω ₁ values of 25, 50, 100, 150, and 200 Hz correspond to B ₁ values of 0.59, 1.17, 2.35, 3.52, and 4.70 μT, respectively. Figure 4(c) shows Z-spectra as a function of Δωoff for various M0b/M0avalues (1/500, 1/250, 1/125, 1/100, and 1/50).

In the above simulations, we assumed that T1aand T2awere 3 s and 100 ms, respectively, and T1b=1sand T2b=15ms[16]. The chemical shift of protons in pool B was set to be 4 ppm. It should be noted that the chemical shift of 4 ppm corresponds to Δωoff of 1192.8 Hz for the magnetic field strength of 7 T. Unless otherwise indicated, kab  + kba was assumed to be 100 Hz. M0aand M0bwere assumed to be 1 and 1/250, respectively. The saturation time and ω ₁ were taken as 2 s and 100 Hz, respectively. The matrix exponential and Kronecker tensor product were calculated using the MATLAB® functions “expm” and “kron,” respectively.

The peaks at 0 Hz (0 ppm) and 1192.8 Hz (4 ppm) in Figure 4 derived from pool A and pool B, respectively. As shown in Figure 4(a) and Figure 4(b) , Z-spectra changed largely depending on the saturation time and ω ₁, that is, Z-spectra became broad and tended to saturate with increasing saturation time and ω ₁. As shown in Figure 4(c) , the peaks at 1192.8 Hz increased with increasing M0b/M0avalue.

Figure 5 shows cases for the three-pool chemical exchange model ( Figure 3 ) consisting of bulk water (pool A) and two labile proton pools (pool B and pool C). In these cases, we assumed that T1a=3s, T2a=100ms, T1b=T1c=1s, and T2b=T2c=15ms[16]. The chemical shifts of two labile proton pools were set to be 4 ppm (Δωoff  = 1192.8 Hz for the magnetic field strength of 7 T) and 5 ppm (Δωoff  = 1491.0 Hz for 7 T). Unless otherwise indicated, kab  + kba , kac  + kca , and kbc  + kcb were assumed to be 100 Hz, 300 Hz, and 100 Hz, respectively. M0a, M0b, and M0cwere assumed to be 1, 1/250, and 1/500, respectively. The saturation time and ω ₁ were taken as 5 s and 50 Hz, respectively.

Figure 5(a) shows Z-spectra as a function of Δωoff for various saturation times (0.5, 1, 2, 5, and 10 s). The peaks at 0 Hz (0 ppm), 1192.8 Hz (4 ppm), and 1491.0 Hz (5 ppm) derive from pool A, pool B, and pool C, respectively. As shown in Figure 5(a) , Z-spectra changed largely depending on the saturation time, that is, Z-spectra became broad and tended to saturate with increasing saturation time. Figure 5(b) shows Z-spectra as a function of Δωoff for various ω ₁ values (25, 50, 100, 150, and 200 Hz). As in Figure 4(b) , Z-spectra became broad with increasing ω ₁ value. Figure 5(c) shows Z-spectra as a function of Δωoff for various M0c/M0avalues (1/500, 1/250, 1/125, 1/100, and 1/50). The peaks at 1491.0 Hz increased with increasing M0c/M0avalue.

Figure 6(a) shows the MTR_asym values given by Eq. (19) as a function of ω ₁ for various saturation times (0.5, 1, 2, 5, and 10 s) in the two-pool chemical exchange model ( Figure 1 ), whereas Figure 6(b) shows those as a function of saturation time for various ω ₁ values (25, 50, 100, 150, and 200 Hz). As shown in Figure 6(a) , when ω ₁ was small, MTR_asym tended to increase with increasing ω ₁ and saturation time. However, when ω ₁ was large, MTR_asym tended to saturate or decrease with increasing ω ₁ value, depending on the saturation time. As shown in Figure 6(b) , MTR_asym tended to saturate with increasing saturation time for all ω ₁ values.

Figure 7(a) shows the PTR values given by Eq. (20) as a function of ω ₁ for various saturation times (0.5, 1, 2, 5, and 10 s) in the two-pool chemical exchange model ( Figure 1 ), whereas Figure 7(b) shows those as a function of saturation time for various ω ₁ values (25, 50, 100, 150, and 200 Hz). As shown in Figure 7 , although PTR showed almost the same tendency with MTR_asym ( Figure 6 ), the change in the PTR value depending on the saturation time or ω ₁ was larger than that in the MTR_asym value.

In this study, we presented a simple equation for solving the time-dependent Bloch-McConnell equations, in which our previous method [18] and the approach presented by Koss et al. [19] were combined. Our method can be easily expanded to multi-pool chemical exchange models by modifying the matrix A in Eq. (2). We previously reported that the solutions obtained by our method agreed with the analytical solutions given by Mulkern and Williams, [21] and the numerical solutions obtained using a fourth/fifth-order Runge–Kutta-Fehlberg (RKF) algorithm [18], indicating the validity of our method. In addition, our method considerably reduced the computation time as compared with the RKF algorithm [18]. These results suggest that our method will be useful in calculating the parameters such as the exchange rate of CEST agents using the non-linear least-squares fitting method [17].

As previously described, the so-called Z-spectrum has usually been used to analyze the CEST effect [18]. The Z-spectrum is obtained by plotting the z magnetization component of bulk water protons (Mza) in the form of Mzaversus Δωoff [Eq. (18)]. Figure 4(a) and Figure 5(a) showed that the saturation time affected the Z-spectra, and the CEST effect increased and saturated with increasing saturation time. The fact that the CEST effect saturates with increasing saturation time is more clearly confirmed by the relationship between MTR_asym or PTR, and the saturation time shown in Figure 6(b) or Figure 7(b) . As shown in Figure 4(b) and Figure 5(b) , ω ₁ also affected the Z-spectra. Although the CEST effect increased with increasing ω ₁ value, the separation among peaks in the Z-spectrum plots degraded with increasing ω ₁ value. The influence of ω ₁ on the CEST effect is also clearly demonstrated by the relationship between MTR_asym and PTR, and ω ₁ shown in Figure 6(a) or Figure 7(a) . The use of large ω ₁ may directly saturate bulk water protons, causing the so-called spillover effect [18]. The results shown in Figures 4 – 7 suggest that the values of imaging parameters in CEST MRI such as the saturation time and ω ₁ must be determined in consideration of both the CEST effect and spillover effect. Our method is useful for determining the optimal values of imaging parameters in CEST MRI.

2.4. Calculation of R _1ρ and R _2ρ

The longitudinal relaxation rate in the rotating frame (R _1ρ ) can be obtained from the negative of the largest (least negative) real eigenvalue (λ ₁) of the matrix A in Eq. (2), that is, R _1ρ  =  − λ ₁ [19, 22].

The transverse relaxation rate in the rotating frame (R _2ρ ) can be obtained from the absolute value of the largest real part of the complex eigenvalue (λ ₂) of the matrix A in Eq. (2), that is, R _2ρ  = |Re(λ ₂)| [22], where Re denotes the real part of a complex number.

Figure 8 shows the common logarithm of R _1ρ (a) and R _2ρ (b) as a function of Δωoff for saturation times of 0.5, 1, 2, 5, and 10 s in the two-pool chemical exchange model ( Figure 1 ). The peaks at 0 Hz (0 ppm) and 1192.8 Hz (4 ppm) derive from pool A and pool B, respectively. As shown in Figure 8 , R _1ρ and R _2ρ were not affected by the saturation time.

Figure 9 shows the common logarithm of R _1ρ (a) and R _2ρ (b) as a function of Δωoff for ω ₁ values of 25, 50, 100, 150, and 200 Hz in the two-pool chemical exchange model ( Figure 1 ). As shown in Figure 9 , both parameters became broad with increasing ω ₁ value.

As described above, R _1ρ and R _2ρ can be obtained from the negative of the largest (least negative) real eigenvalue and the absolute value of the largest real part of the complex eigenvalue of the matrix A in Eq. (2), respectively. We previously reported that there was good agreement between the R _1ρ and R _2ρ values thus obtained and those obtained numerically [22]. These results appear to indicate the validity of these procedures.

As shown in Figure 8 , R _1ρ and R _2ρ were not affected by the saturation time, because the matrix A in Eq. (2) is independent of the saturation time. When ω ₁ was varied, the influence of ω ₁ on R _1ρ and R _2ρ increased with increasing ω ₁ value ( Figure 9 ). Especially, the separation between peaks in the R _1ρ plots degraded with increasing ω ₁ value [ Figure 9(a) ]. This also appears to be due to the spillover effect.

3.1. Principle of spin-locking

Longitudinal relaxation time in the rotating frame (T _1ρ ) has been demonstrated to be effective for probing the slow-motion interactions between motion-restricted water molecules and their local macromolecular environment [23] and provides novel image contrast that is not available from conventional MRI techniques. The imaging of biologic tissue based on T _1ρ is currently being investigated for various tissues, including articular cartilage, breast, and head and neck [24, 25, 26]. In T _1ρ -weighted MRI of tissues, the image is sensitive to molecular processes that occur over a range of frequencies determined by the amplitude of an applied SL pulse [23].

As pointed out by Jin et al. [27], the SL approach is useful for improving the signal-to-noise ratio (SNR) in CEST MRI. Furthermore, Kogan et al. [28] demonstrated that a combination of the CEST and SL approaches is useful for detecting proton exchange in the slow-to intermediate-exchange regimes.

As earlier described, the Bloch-McConnell equations for the two-pool chemical exchange model ( Figure 1 ) in the rotating frame with the same frequency as that of the RF-pulse irradiation is given by Eq. (2) [18, 29]. The solution of Eq. (2) can be given by [18]

Figure 10 illustrates the image of the pulse sequence with SL. We assume that the SL pulse (frequency: ω, amplitude: ω ₁, and frequency offset: Ω) is applied on the x-axis ( Figure 2 ). The effective magnetic field (B1eff)anditsangle with respect to thez−axis(θ) are given by B1eff=ω12+Ω2/γand θ = tan⁻¹(ω ₁/Ω), respectively ( Figure 2 ). To achieve SL, the magnetization is first flipped by the θ-degree RF pulse (frequency: ω and amplitude: ω1θ) to the x-z plane, then spin locked by B1effforaduration oftSL, and then flipped back to the z-axis for imaging ( Figure 10 ). The θ-degree RF pulse for flipping is applied on the –y axis, that is, ϕ =  − π/2, whereas the θ-degree RF pulse for flipping back is applied on the y axis, that is, ϕ = π/2. The θ-degree rotation matrix for flipping [R(θ)] is given by [30].

where ω1θand tθ denote the amplitude and the duration of the θ-degree RF-pulse irradiation, respectively ( Figure 10 ), and ω1θ×tθ=θ. Thus, we obtain the magnetization vector immediately after SL for a duration of tSL [M ⁻(tSL )] as.

The θ-degree rotation matrix for flipping back to the z-axis [R(−θ)] is given by.

Thus, the magnetization vector after flipping back to the z-axis [M ⁺(tSL )] is given by.

Note that Ω and θ are taken to be 0 and π/2, respectively, for an on-resonance SL sequence, whereas the saturation pulse is applied without flipping the magnetization in the sequence without SL such as the conventional CEST sequence [15]. Therefore, the magnetization vector after the saturation pulse [M(tSAT )] in the conventional CEST MRI is simply expressed as.

where tSAT denotes the duration of saturation.

3.2. Calculation of T _1ρ

T _1ρ can be obtained numerically by fitting the z component of magnetization for tSL [M ⁺(tSL ) given by Eq. (25)] in pool A [MzatSL] to the following equation [30]:

where Mzssadenotes the steady-state z component of magnetization in pool A. In this study, we used the Simplex method [31] to calculate T _1ρ from Eq. (27).

The approximate solution for T _1ρ has been derived by Trott and Palmer [29]:

where θ = tan⁻¹(ω ₁/Ω), Rex=PaPbΔω2kex/ωae2ωbe2/ωe2+kex2, ωae=ω12+Δωa2, ωbe=ω12+Δωb2, ωe=ω12+Ω2, Ω=ω¯−ω, ω¯=Paωa+Pbωb, Δω = Δωb  − Δωa  = ωb  − ωa , and kex  = kab  + kba . Pa and Pb are the fractional sizes of pool A and pool B, and are given by Pa=M0a/M0a+M0band Pb=M0b/M0a+M0b, respectively. R ₁ and R ₂ are the population-averaged relaxation rates, and are given by R1=PaR1a+PbR1band R2=PaR2a+PbR2b, respectively. It should be noted that Ω is the population-averaged offset frequency in this case. Mzssain Eq. (27) is approximated by [27].

Figure 11 shows an example of the three-dimensional plots of the magnetization vector in pool A in the two-pool chemical exchange model ( Figure 1 ). Figure 11(a) and 11(b) show cases without and with SL, respectively. In these cases, the relaxation time constants were assumed to be T1a=1.5s, T2a=0.06s, T1b=0.77s, and T2b=0.033s[32]. tθ in Eq. (22) and (24) was taken as 200 μs [27]. ω1θin Eq. (22) and (24) was calculated from ω1θ=θ/tθ. Unless specifically stated, Δω (= ωb  − ωa ) and ω ₁ were assumed to be 2400 and 1000 Hz, respectively. Ω was assumed to be 2000 Hz. Thus, θ was tan⁻¹(ω ₁/Ω) = tan⁻¹(1000/2000) ≈ 26.6 degrees. kex (= kab  + kba ) was assumed to be 1500 Hz, and kab was assumed to be given by kab=M0b/M0akba[18]. M0b/M0awas assumed to be 0.03. As shown in Figure 11(a) , when the SL pulse was not applied, the magnetization vector rotated largely around B1eff. On the other hand, when the SL pulse was applied [ Figure 11(b) ], the magnetization vector moved along B1eff, and the rotation around B1effwas suppressed.

When M0b/M0awas 0.003, there was good agreement between the T _1ρ values calculated from Eq. (27) and Trott and Palmer’s solutions given by Eq. (28) (data not shown). When M0b/M0awas 0.03, some difference was observed between them in the off-resonance case. When M0b/M0awas 0.3, large differences were observed between them in both the on- and off-resonance cases [30].

In this study, we developed a simple and fast method for calculating the magnetization vector in SL CEST MRI, in which a simple matrix equation was derived for solving the time-dependent Bloch-McConnell equations in SL MRI [Eq. (25)] and the θ-degree rotation matrix [Eq. (22)] was introduced for considering the effect of the θ-degree RF pulse for flipping the magnetization. As shown in Figure 11 , the trajectory of the magnetization vector in the sequence with SL could be visualized by calculating M ⁻(tSL ) using Eq. (23), whereas that in the sequence without SL could be visualized by calculating M(tSAT ) using Eq. (26). Although Figure 11 shows the three-dimensional plots observed from one direction, we can observe the trajectory of the magnetization vector from various directions by rotating the plot. If we compared the three-dimensional plots with and without SL ( Figure 11 ), then the effect of SL is well understood. Therefore, our method is helpful for visually understanding the effect of SL. In addition, as our method allows us to simply and quickly calculate the time evolution of the magnetization vector under various study conditions in SL CEST MRI, our method can also be useful for optimizing the study conditions in SL CEST MRI.

As previously described, when M0b/M0awas small, that is, when the population of two pools was highly asymmetric, the T _1ρ values calculated from Eq. (27) agreed with the solutions given by Eq. (28). However, the difference between them increased with increasing M0b/M0a[30]. This finding appears to be due to the fact that Trott and Palmer’s solution [Eq. (28)] was derived by approximating the parameters such as relaxation rates using their population-averaged values, and thus the validity of this approximation decreases with decreasing asymmetry in the populations of the two pools.

Although we treated the two-pool chemical exchange model ( Figure 1 ) for analyzing T _1ρ or R _1ρ in SL CEST MRI, recent investigations have shown the importance of improved theoretical approaches for describing multi-site chemical exchange phenomena [33, 34]. Thus, Trott and Palmer [33] have tried to generalize their approach for T _1ρ or R _1ρ [29]. For such purposes, it is necessary to expand the Bloch-McConnell equations to those based on multi-pool chemical exchange models. Our method can be easily expanded to multi-pool chemical exchange models by modifying the matrix A given by Eq. (4) [20] as previously described, and it is helpful for testing the validity of newly developed approaches for analyzing multi-site chemical exchange phenomena.