Abstract
Chemical exchange saturation transfer (CEST) is one of the contrast mechanisms in magnetic resonance imaging (MRI) and has been used to detect dilute proteins through the interaction between bulk water and labile solute protons. Amide proton transfer (APT) MRI has been developed for imaging diseases such as acute stroke. Moreover, various CEST agents have been explored to enhance the CEST effect. The contrast mechanism of CEST or APT MRI, however, is complex and depends not only on the concentration of amide protons or CEST agents and exchange properties, but also varies with imaging parameters such as radiofrequency (RF) power and magnetic field strength. When there are multiple exchangeable pools within a single CEST system, the contrast mechanism of CEST becomes even more complex. Numerical simulations are useful and effective for analyzing the complex contrast mechanism of CEST and for investigating the optimal imaging parameter values. In this chapter, we present the basics of CEST or APT MRI and a simple and fast numerical method for solving the time-dependent Bloch-McConnell equations for analyzing the behavior of magnetization and/or contrast mechanism in CEST or APT MRI. We also present a method for analyzing the behavior of magnetization in spin-locking CEST MRI.
Keywords
- Bloch-McConnell equations
- numerical solution
- chemical exchange saturation transfer (CEST) MRI
- amide proton transfer (APT) MRI
- spin-locking
1. Introduction
Chemical exchange saturation transfer (CEST) is one of the contrast mechanisms in magnetic resonance imaging (MRI) [1] and has been increasingly used to detect dilute proteins through the interaction between bulk water protons and labile solute protons [2, 3, 4]. Amide proton transfer (APT) MRI has been developed for imaging diseases such as acute stroke and cancer, and is now under intensive evaluation for clinical translation [5, 6]. APT MRI is a particular type of CEST MRI that specifically probes labile amide protons of endogenous mobile proteins and peptides in tissue [5, 6]. In addition to APT MRI [5, 6], useful CEST MRI contrast for clinical imaging can be generated from amine protons [7], hydroxyl protons [8], glycosaminoglycans [9], and glutamate [10], as well as from changes in creatine and lactate concentrations [11]. Glucose and iopamidol have been used as exogenous CEST agents that have been administered to patients [12, 13]. Moreover, various CEST agents have been energetically developed to detect the parameters that reflect tissue molecular environment such as hydrogen ion exponent (pH) and/or to enhance the CEST effect [14].
In CEST or APT MRI, the exchangeable proton spins are saturated, and the saturation is transferred upon chemical exchange to the bulk water pool [1, 15]. As a result, a large contrast enhancement in bulk water can be achieved. The contrast mechanism of CEST or APT MRI, however, is complex and depends not only on the concentration of amide protons or CEST agents, relaxation, and exchange properties but also varies with imaging parameters such as radiofrequency (RF) power and magnetic field strength [15]. When there are multiple exchangeable pools within a single CEST system, the contrast mechanism of CEST becomes all the more complex [16]. Numerical simulations are useful and effective for analyzing the complex CEST contrast mechanism and for investigating the optimal imaging parameter values [17, 18]. In order to perform extensive numerical simulations for CEST or APT MRI, it requires the development of a simple and fast numerical method for obtaining the solutions to the time-dependent Bloch-McConnell equations.
In this chapter, we present the basics of CEST or APT MRI and a simple and fast numerical method for solving the time-dependent Bloch-McConnell equations for analyzing the behavior of magnetization and/or contrast mechanism in CEST or APT MRI. We also present it in SL CEST MRI.
2. Bloch-McConnell equations in the presence of CEST
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].

Figure 1.
Two-pool chemical exchange model.
where superscripts a and b show the parameters in pool A and pool B, respectively. For example,

Figure 2.
Parameters for analyzing the behavior of magnetization in the rotating frame.
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( ω , ω 1, 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( ω , ω 1, 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].

Figure 3.
Three pool chemical exchange model.
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, MTRasym, 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
The magnetization transfer asymmetry (MTRasym) analysis has been performed using the following equation [18]:
Instead of MTRasym, the following equation for proton transfer ratio (PTR) has also been used for analyzing the CEST effect [18]:
where
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 ω 1 values (25, 50, 100, 150, and 200 Hz). It should be noted that because B 1 = ω 1/γ, ω 1 values of 25, 50, 100, 150, and 200 Hz correspond to B 1 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

Figure 4.
Z-spectra as a function of Δ ω off for various values of saturation time a , ω 1 b , and M 0 b / M 0 a (c) in the two-pool chemical exchange model.
In the above simulations, we assumed that
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 ω 1, that is, Z-spectra became broad and tended to saturate with increasing saturation time and ω 1. As shown in Figure 4(c) , the peaks at 1192.8 Hz increased with increasing
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

Figure 5.
Z-spectra as a function of Δ ω off for various values of saturation time a , ω 1 b , and M 0 c / M 0 a (c) in the three-pool chemical exchange model.
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 ω 1 values (25, 50, 100, 150, and 200 Hz). As in Figure 4(b) , Z-spectra became broad with increasing ω 1 value. Figure 5(c) shows Z-spectra as a function of Δωoff for various
Figure 6(a) shows the MTRasym values given by Eq. (19) as a function of ω 1 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 ω 1 values (25, 50, 100, 150, and 200 Hz). As shown in Figure 6(a) , when ω 1 was small, MTRasym tended to increase with increasing ω 1 and saturation time. However, when ω 1 was large, MTRasym tended to saturate or decrease with increasing ω 1 value, depending on the saturation time. As shown in Figure 6(b) , MTRasym tended to saturate with increasing saturation time for all ω 1 values.

Figure 6.
(a) MTRasym values as a function of ω 1 for various saturation times in the two − pool chemical exchange model. (b) MTRasym values as a function of saturation time for various ω 1 values.
Figure 7(a) shows the PTR values given by Eq. (20) as a function of ω 1 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 ω 1 values (25, 50, 100, 150, and 200 Hz). As shown in Figure 7 , although PTR showed almost the same tendency with MTRasym ( Figure 6 ), the change in the PTR value depending on the saturation time or ω 1 was larger than that in the MTRasym value.

Figure 7.
(a) PTR values as a function of ω 1 for various saturation times in the two − pool chemical exchange model. (b) PTR values as a function of saturation time for various ω 1 values.
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 (
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 (λ 1) of the matrix A in Eq. (2), that is, R 1ρ = − λ 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 (λ 2) of the matrix A in Eq. (2), that is, R 2ρ = |Re(λ 2)| [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 8.
(a) Common logarithm of R 1 ρ and (b) R 2 ρ values as a function of Δω off for various saturation times in the two-pool chemical exchange model.
Figure 9 shows the common logarithm of R 1ρ (a) and R 2ρ (b) as a function of Δωoff for ω 1 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 ω 1 value.

Figure 9.
(a) Common logarithm of R 1 ρ and (b) R 2 ρ values as a function of Δω off for various ω 1 values in the two-pool chemical exchange model.
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 ω 1 was varied, the influence of ω 1 on R 1ρ and R 2ρ increased with increasing ω 1 value ( Figure 9 ). Especially, the separation between peaks in the R 1ρ plots degraded with increasing ω 1 value [ Figure 9(a) ]. This also appears to be due to the spillover effect.
3. Spin-locking CEST MRI
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: ω 1, and frequency offset: Ω) is applied on the x-axis ( Figure 2 ). The effective magnetic field (

Figure 10.
Diagram of spin-locking pulse sequence.
where
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 [
where
The approximate solution for T 1ρ has been derived by Trott and Palmer [29]:
where θ = tan−1(ω 1/Ω),
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

Figure 11.
Three-dimensional plots of the magnetization vector in pool A in the two-pool chemical exchange model. (a) and (b) show cases without and with spin-locking pulse, respectively.
When
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
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.
4. Correction of B0 and B1
As previously described, the CEST effect has usually been analyzed using MTRasym [Eq. (19)] or PTR [Eq. (20)]. However, these parameters are susceptible to the B0 inhomogeneity of the static magnetic field. When there exists the B0 inhomogeneity, the spillover effect is no longer symmetric. Furthermore, the B1 inhomogeneity of the RF pulse may also cause spatial variation in labeling efficiency and spillover factor [35]. Apart from the efforts in improving magnetic field inhomogeneities using hardware-based methods, such as parallel transmit technologies [36], post-processing algorithms have been developed for field inhomogeneity correction [37, 38].
Kim et al. [37] showed that direct water saturation imaging allows measurement of the absolute water frequency in each voxel, allowing proper centering of Z-spectra on a voxel-by-voxel basis independent of spatial B 0 field variations, and that the B0 inhomogeneity in CEST MRI can be corrected on a voxel-by-voxel basis through the centering of Z-spectra. This method is called “water saturation shift referencing (WASSR)” approach. This method, however, would require acquisition of saturation images at 20–40 frequencies [38]. Since the SNR of CEST MRI is low, multiple acquisitions for each frequency offset of complete Z-spectra would be needed, which is not practical in the clinical setting. Zhou et al. demonstrated that a practical six-offset multi-acquisition method combined with a single reference Z-spectrum to acquire high-SNR CEST MRI can accomplish improved CEST MRI with B0 inhomogeneity correction within an acceptable scanning time [38].
A B1-correction of CEST contrasts is crucial for the evaluation of data obtained in clinical studies at high field strengths with strong B1-inhomogeneities. To correct for the B1 inhomogeneity, a B1 map is acquired for correction of Z-spectra using either a calibration [39] or an interpolation approach [40]. Singh et al. [39] developed an approach for B 1 inhomogeneity correction based on acquiring calibration data at a coarsely sampled B 1 values in conjunction with the measured B 1 maps, whereas Windschuh et al. [40] developed an approach based on Lorentzian line fits.
The comprehensive methods like simultaneous mapping of B0 and B1 fields [35, 41], and model-based correction algorithm, [42] have also been developed to improve the accuracy of MTRasym or PTR.
Acknowledgments
This work was supported, in part, by a Grant-in-Aid for Challenging Exploratory Research (Grant No. 25670532) from the Japan Society for the Promotion of Science.