Numerical Analysis of the Magnetization Behavior in Magnetic Resonance Imaging in the Presence of Multiple Chemical Exchange Pools Kenya

The purpose of this study was to demonstrate a simple and fast method for solving the time-dependent Bloch-McConnell equations describing the behavior of magnetization in magnetic resonance imaging (MRI) in the presence of multiple chemical exchange pools. First, the time-dependent BlochMcConnell equations were reduced to a homogeneous linear differential equation, and then a simple equation was derived to solve it using a matrix operation and Kronecker tensor product. From these solutions, the longitudinal relaxation rate (R1ρ) and transverse relaxation rate in the rotating frame (R2ρ) and Z-spectra were obtained. As illustrative examples, the numerical solutions for linear and star-type three-pool chemical exchange models and linear, startype, and kite-type four-pool chemical exchange models were presented. The effects of saturation time (ST) and radiofrequency irradiation power (ω1) on the chemical exchange saturation transfer (CEST) effect in these models were also investigated. Although R1ρ and R2ρ were not affected by the ST, the CEST effect observed in the Z-spectra increased and saturated with increasing ST. When ω1 was varied, the CEST effect increased with increasing ω1 in R1ρ, R2ρ, and Z-spectra. When ω1 was large, however, the spillover effect due to the direct saturation of bulk water protons also increased, suggesting that these parameters must be determined in consideration of both the CEST and spillover effects. Our method will be useful for analyzing the complex CEST contrast mechanism and for investigating the optimal conditions for CEST MRI in the presence of multiple chemical exchange pools.


Introduction
Chemical exchange saturation transfer (CEST) is a novel contrast mechanism for magnetic resonance imaging (MRI) [1] and has been increasingly used to detect dilute proteins via the interaction between labile solute protons and bulk water protons [2] [3] [4].Moreover, amide proton transfer (APT) imaging, a particular type of CEST MRI that specifically probes labile amide protons of endogenous mobile proteins and peptides in tissue, has been explored for imaging diseases such as acute stroke and tumor, and is currently under intensive evaluation for clinical translation [5] [6].Furthermore, various CEST agents have been actively developed to detect the parameters that reflect tissue pH and molecular environment and/or to enhance the CEST effect [7].However, CEST or APT MRI contrast mechanism is complex, depending on not only the concentration of CEST agents or amide protons, exchange and relaxation properties, but also varying with experimental conditions such as magnetic field strength and radiofrequency (RF) power [8].When there are multiple exchangeable sites within a single CEST system, the CEST contrast mechanism becomes even more complex [9].Thus, in analyzing the complex CEST contrast mechanism and for investigating the optimal study conditions, numerical simulations are useful and effective [10] [11].To perform extensive numerical simulations for CEST or APT MRI, it will be necessary to develop a simple and fast method for obtaining the numerical solutions to the time-dependent Bloch-McConnell equations.
The purpose of this study was to present a simple and fast method for solving the time-dependent Bloch-McConnell equations for analyzing the behavior of magnetization in MRI in the presence of multiple chemical exchange pools.ω and b ω are the Larmor frequencies in pool A and pool B, respectively, and ω is the frequency of the RF irradiation applied along the x axis of the rotating frame.ω 1 is the amplitude of the RF irradiation.

Bloch-McConnell Equations in a Two-Pool Chemical Exchange Model
The differential equations given by Equation ( 1) can be combined into one vector equation (homogeneous linear differential equation) [11]: where where T in Equation (3) denotes the matrix transpose.According to Koss et al. [12], the matrix A can be given by 0 0 where E is the evolution matrix [12] and C is the constant-term matrix.Furthermore, E is given by In the case of A given by Equation ( 4), R is reduced to where and K in Equation ( 6) is given by where I is a 3-by-3 identity matrix and ⊗ denotes the Kronecker tensor product.C in Equation ( 5) is given by [ ] The solution of Equation ( 2) can be given by [11] ( ) ( ) where t represents the so-called saturation time and ( ) is the matrix of initial values at 0 t = .

Linear Three-Pool Chemical Exchange Model
Figure 2(a) illustrates a linear three-pool chemical exchange model in which pool A represents the bulk water pool.In this case, R and K are given by [12] and 0 0 respectively.c R in Equation ( 13) is given by Equation ( 8) in which the sub- script a and superscript a are replaced by c. C is given by pool A represents the bulk water pool.In this case, K is given by [12] ab ac ba ca A represents the bulk water pool.In this case, R and K are given by [12]

Linear Four-Pool Chemical Exchange Model
and respectively.R d in Equation ( 17) is given by Equation (8) in which the subscript a and superscript a are replaced by d.C is given by

Kite-Type Four-Pool Chemical Exchange Model
It should be noted that mass balance imposes the following relationship be- ( ) and ( )

Calculation of R1ρ, R2ρ, and Z-Spectra
The longitudinal relaxation rate in the rotating frame ( ) R ρ was obtained from the negative of the largest (least negative) real eigenvalue ( ) The transverse relaxation rate in the rotating frame ( ) 2 R ρ was obtained from the absolute value of the largest real part of the complex eigenvalue ( ) Re R ρ λ = [13].
The CEST effect has usually been analyzed using the so-called Z-spectrum [11].Thus, we calculated the Z-spectrum by using the following equation: where ( )

Simulation Studies
Because we have already treated a two-pool chemical exchange model in our previous paper [11], we treated three-pool and four-pool chemical exchange models in this study.
First, we considered the three-pool exchange model consisting of bulk water (pool A) and two labile proton pools (pool B and pool C) as illustrative examples.
In this case, we assumed that the longitudinal ( ) For four-pool exchange models, we simulated one nuclear Overhauser effect site (pool D) in addition to the above bulk water (pool A) and two labile proton pools (pool B and pool C).We assumed that the longitudinal ( )

Results
Figure 4 shows the 1 R ρ (a), 2 R ρ (b), and Z-spectra (c) as a function of offset frequency ( ) off ω ∆ for various saturation times (0.5, 1, 2, 5, and 10 s) in the linear three-pool chemical exchange model (Figure 2(a)).It should be noted that the common logarithm of the 1 R ρ value was plotted in Figure 4(a) in order to enlarge the change in the 1 R ρ value.The peaks at 0 Hz (0 ppm), 1192.8Hz (4 ppm), and 1491.0Hz (5 ppm) derive from pool A, pool B, and pool C, respectively.As shown in Figure 4, 1 R ρ and 2 R ρ were not affected by the saturation time, whereas Z-spectra changed largely depending on the saturation time, i.e., Z-spectra became broad and tended to saturate with increasing saturation time.
Figure 5 shows the 1 R ρ (a), 2 R ρ (b), and Z-spectra (c) as a function of

Discussion
In this study, we developed a simple equation for solving the time-dependent Bloch-McConnell equations in the presence of multiple chemical exchange pools by combining our previous method [11] and the approach presented by Koss et al. [12].As described in our previous paper [11], the numerical solutions ob- tained by our method agreed with the analytical solutions given by Mulkern and Williams [14].We also compared the solutions obtained by our method for the two-pool exchange model with those obtained using a fourth/fifth-order Runge-Kutta-Fehlberg (RKF) algorithm and found that there was a good agreement between them [11].These results appear to indicate the validity of our method.
Furthermore, the computation time was considerably reduced when using our method (by a factor of approximately 2500 compared to the case when using the RKF algorithm [11]).Thus, our method can be included in the nonlinear leastsquares fitting routine to calculate parameters such as the exchange rate or lifetime of CEST agents [10].
For calculating the solutions to the time-dependent Bloch-McConnell equations using Equation ( 12), most computation time is spent calculating the eigenvectors and eigenvalues of matrix A. However, it is necessary to carry out this calculation only once regardless of t in Equation (12).As previously described, in this study, the matrix exponential was computed using the MATLAB ® function "expm", in which a scaling and squaring algorithm with Pade approximation [15] has been used.
In our previous study [11], we used the two-pool exchange model for CEST or APT MRI as an illustrative example.As pointed out by Woessner et al. [10], paramagnetic CEST agents often have more than one type of exchangeable proton.
For such cases, it is necessary to expand the Bloch-McConnell equations to multi-pool exchange models.Recently, Koss et al. [12] presented a generalized expression for the evolution matrix in the Bloch-McConnell equations in the presence of multiple chemical exchange sites.In their method, the Kronecker tensor K. Murase product was used [12].As shown in this study, our method could be easily extended to multi-pool chemical exchange models by modifying matrix A in Equation (12) with use of their approach.
As previously described, 1 R ρ was obtained from the negative of the largest (least negative) real eigenvalue of matrix A in Equation (2). 2 R ρ was obtained from the absolute value of the largest real part of the complex eigenvalue of matrix A in Equation ( 2).We previously compared the 1 R ρ and 2 R ρ values thus obtained with those obtained numerically and found that there was a good agreement between them [13], indicating the validity of these procedures.
The spectral dependence of CEST is determined by sweeping the RF irradiation frequency while monitoring the water resonance [1].As previously described, the CEST effect has usually been analyzed using the so-called Z-spectrum [11].The Z-spectrum is obtained by plotting the z component of the magnetization of pool A, i.e., bulk water proton ( )  2) is independent of the saturation time.On the other hand, the Z-spectra were affected by the saturation time, i.e., the CEST effect observed in the Z-spectra increased and saturated with increasing saturation time (Figure 4(c)).When the RF irradiation power ( ) ω was varied, the CEST effect increased with increasing 1 ω in 1 R ρ , 2 R ρ , and Z-spectra (Figure 5 and Figure 7).When 1 ω is large, however, its saturation bandwidth is broad and thus may directly saturate the bulk water, causing the so-called spillover effect [8].As shown in Figure 5 and Figure 7, the separation among peaks in the 1 R ρ and Z-spectrum plots degraded with increasing 1 ω , which appears to be due to the spillover effect.These results suggest that 1 ω must be determined in con- sideration of both the CEST effect and the spillover effect.Simulation studies with use of our method will be useful especially in such a case.
As shown in Figure 6 R value in a two-pool chemical exchange model [13].This would also be applicable to the case of multi-pool chemical exchange models.Thus, the above finding would be able to be explained by the fact that the Recently, Koss et al. [12] presented analytical expressions for 1 R ρ in the presence of multiple chemical exchange sites and pointed out that analytical solutions facilitate understanding of the relationship between model parameters and the phenomenological relaxation rate constant and can lead to new methodological advances.Numerical solutions also appear to be useful for optimizing parameters such as the saturation time and 1 ω for acquiring CEST or APT MRI data, because they can be easily obtained under various and/or complex study conditions in which analytical solutions may not always be obtained.

Conclusion
We presented a simple and fast numerical method for solving the time-dependent Bloch-McConnell equations in the presence of multiple chemical exchange pools by combining our previous method [11] and the approach presented by Koss et al. [12].The present method will be useful for analyzing the complex CEST con- trast mechanism and for investigating the optimal conditions for CEST MRI in the presence of multiple chemical exchange pools.

Figure 1
Figure 1 illustrates a two-pool chemical exchange model in which pool A represents the bulk water pool.The time-dependent Bloch-McConnell equations in the two-pool exchange model for CEST MRI are given by[10] [11]

Figure 1 . 1 a R and 2 aRM
Figure 1.Illustration of a two-pool chemical exchange model.k ab and k ba represent the exchange rates from pool A to pool B and from pool B to pool A, respectively.

teA
is the matrix exponential.

Figure 2 .
Figure 2. Illustration of three-pool chemical exchange models.(a) and (b) show linear and star-type three-pool chemical exchange models, respectively.As in the case of k ab , k ac , k ca , k bc , and k cb represent the exchange rates from pool A to pool C, from pool C to pool A, from pool B to pool C, and from pool C to pool B, respectively.

Figure 2 (
Figure 2(b) illustrates a triangular three-pool chemical exchange model in which

Figure 3 (
Figure 3(a) illustrates a linear four-pool chemical exchange model in which pool

Figure 3 (
Figure 3(b) illustrates a star-type four-pool chemical exchange model in whichpool A represents the bulk water pool.In this case, K is given by[12].

Figure 3 .
Figure 3. Illustration of four-pool chemical exchange models.(a), (b), and (c) show linear, star-type, and kite-type four-pool chemical exchange models, respectively.As in the case of k ab , k ad , k da , k cd , and k dc represent the exchange rates from pool A to pool D, from pool D to pool A, from pool C to pool D, and from pool D to pool C, respectively.

Figure 3 (
Figure 3(c) illustrates a kite-type four-pool chemical exchange model in which pool A represents the bulk water pool.Although there are no chemical exchanges between pool C and pool D in the star-type four-pool chemical exchange model (Figure 3(b)), there are exchanges between them in the kite-type model (Figure3(c)).In this case, K is given by[12] tween the exchange rates ( ij k and ) ji k of pool I and pool J [10]: and 100 ms, respectively, and were 1 s and 15 ms for two labile protons, i.e., chemical shifts for two labile protons were set to be 4 ppm and 5 ppm.It should be noted that the chemical shifts of 4 ppm and 5 ppm correspond to off ω ∆ of 1192.8Hz and 1491.0Hz, respectively, for the magnetic field strength of 7 T.Unless specifically stated, to be 1, 1/250, and 1/500, respectively.The saturation time and 1 ω were taken as 5 s and 50 Hz, respectively.
and 5 ms, respectively [9].The chemical shift for pool D was set to be −3.5 ppm.It should be noted that the chemical shift of −3.5 ppm corresponds to off ω ∆ of −1043.7 Hz for the magnetic field strength of 7 T. Unless specifically stated, be 10 Hz and 10 Hz, respectively.
were performed using MATLAB® (The MathWorks Inc., Natick, MA, USA) on an Intel Core TM i7-4790 CPU (3.6 GHz) with 8-GB RAM.The matrix exponential and Kronecker tensor product were calculated using the MATLAB® functions "expm" and "kron", respectively.

Figure 6 M∆ for various 1 ωFigure 4 .
Figure 6 shows the 1 R ρ (a), 2 R ρ (b), and Z-spectra (c) as a function of off ω ∆ for various

Figure 8 shows the 1 R
Figure 8 shows the 1 R ρ (a), 2 R ρ (b), and Z-spectra (c) as a function of off ω ∆ for various

Figure 9 shows the 1 R 2 R
Figure 9 shows the 1 R ρ (a), 2 R ρ (b), and Z-spectra (c) as a function of off ω ∆ for various cd

Figure 9 .M
Figure 9. (a) 1 R ρ , (b) 2 R ρ , and (c) Z-spectra as a function of offset frequency ( ) off ω ∆ for various cd dc k k + values (10, 100, 200, 500, and 1000 Hz) in the kite-type four-pool chemical exchange model (Figure 3(c)).In these simulations, saturation time, 1 ω , 0 a M , 0 b M , 0 c M , and 0 d M were assumed to be 5 s, 50 Hz, 1, 1/250, 1/500, and 1/500, respectively.ab ba k k + , ac ca k k + , and Figure 4(b), 1R ρ and 2 R ρ were not affected by the saturation time, because the matrix A in Equation (2) is independent of the saturation time.On the other hand, the Z-spectra were affected by the saturation time, i.e., the CEST effect

M
(b) and Figure 8(b), the 2 R ρ values increased with in- increased, the interaction between bulk water protons and protons in pool C or pool D would increase, leading to an increase of 2 a R .Our previous study demonstrated that the 2 R ρ value increased with increasing 2 a