0760

Simulation of Bloch and Bloch-McConnell equations - speed and accuracy
Christina Graf1, Armin Rund2, Christoph Stefan Aigner1, Karl Kunisch2, and Rudolf Stollberger1

1Institute of Medical Engineering, Graz University of Technology, Graz, Austria, 2Institute of Mathematics and Scientific Computing, Karl-Franzens-University Graz, Graz, Austria

Synopsis

Accurate Bloch and Bloch-McConnell model simulation is essential for proper modelling and becomes crucial for iterative RF pulse design methods. This work demonstrates the advantages of a piecewise analytical solution using eigenvalues and eigenvectors and an asymmetric and symmetric operator splitting scheme. Those three solvers are compared to five other numerical solution methods regarding accuracy and runtime. The symmetric operator splitting prevails in both compared to the others. Moreover, it is quadratic convergent with respect to the time step and therefore seems to be a good choice for optimal control approaches for Bloch and for Bloch-McConnell models.

INTRODUCTION

Iterative RF pulse optimization1,2 requires multiple Bloch or extended Bloch model simulations with the need for fast and accurate solution methods. Especially in the case of chemical exchange saturation transfer (CEST) the typically long pulse durations are a challenge for iterative optimization methods3. In this work, different analytical and numerical methods for solving the Bloch and Bloch-McConnell equations are derived and compared in extensive numerical studies regarding accuracy and runtime.

THEORY

Both the Bloch and Bloch-McConnell equations can mathematically be described as a system of differential equations $$\dot{M}(t) = A(t) \cdot M(t) + b, \: \: \: M(0)=M^{0} \: \: \text{for} \: \: t \in (0,T).$$

Therein, M(t) is the magnetization vector, M0 the initial magnetization and b the inhomogeneity including the longitudinal relaxation time. The coefficient matrix A(t) is non-constant, therefore a full analytic solution does not exist in general. However, if a piecewise constant discretization is chosen for RF pulse and slice-selective gradient (Gs) amplitude, an analytical solution per time step can then be derived using eigenvalues and eigenvectors (EV). On the other hand, by splitting the system matrices into relaxation/chemical exchange and rotation parts, an asymmetric (ASY)4 and a symmetric (SY) operator splitting scheme can be defined, where all subsystems are solved exact. The splitting error is known to be of first-order for asymmetric and of second-order for symmetric operator splitting.

METHODS

The three methods (EV, ASY and SY) are implemented in MATLAB (The MathWorks, Inc., Natick, USA) and are analyzed in extensive numerical studies, both for the Bloch equations and the Bloch-McConnell equations. For the former, the methods are compared to five other numerical solution methods including the spin domain solution with neglected relaxation terms SD5, the MATLAB built-in functions ode45 and expM, an implicit Euler scheme IE and the Crank-Nicholson method CN. The numerical comparison is performed for a complex 180o RF pulse and a trapezoidal Gs shape with a duration of T=1.9ms discretized with 1104 points. The spatial domain is set to (-50,50)mm using 101 points. Different relaxation times are investigated, T1=T2=$$$\infty$$$, T1=1331ms, T2=80ms, and T1=400ms, T2=5ms. The proposed Bloch-McConnell simulation methods (ASY and SY) are compared to the analytical EV solver for a series of Gaussian and block RF pulses with a total duration of T=2100ms discretized with 4421 points in the first place for an amide proton system in white matter6. Later, a convergence study with respect to time discretization is performed. The z-spectrum (-5,5)ppm is discretized using 101 points. The relaxation times of the water proton pool are set to T1,w=1048ms and T2,w=69ms and of the solute proton pool to T1,s=1048ms and T2,s=15ms.

RESULTS AND DISCUSSION

Figure 1 and 2 show the accuracy of all Bloch solvers. In absence of relaxation effects, most solvers yield exact results in machine precision, larger errors are observed for CN and IE. A closer look reveals that IE results in a magnetization degradation over time (not shown). Furthermore, we analyze the behavior in the presence of relaxation effects. IE, CN, expM and SD yield larger errors, the latter since it does not account for relaxation effects. In contrast, ode45, ASY and SY perform with higher precision. Regarding runtime, ASY is slightly faster than SY, which in turn is 3 orders of magnitude faster than ode45 (for the chosen accuracy), and EV. In general, SY outperforms ASY in the accuracy for the chosen time-step. This is also true in presence of chemical exchange saturation transfer shown in Figure 3, 4 and 5. Figure 5 lists the numerically verified linear and quadratic convergence with respect to the time step for ASY and SY in high precision for the series of RF block pulses.

CONCLUSION

We presented a piecewise analytical solver EV for the Bloch and Bloch-McConnell equations that allows exact simulation at the cost of an increased runtime. The symmetric operator splitting SY shows an excellent accuracy together with quadratic convergence in the time-step. Since its runtime is only slightly increased compared to ASY, it seems to be a good choice for optimal control approaches for Bloch as well as for Bloch-McConnell models.

Acknowledgements

Supported by the Austrian Science Fund (FWF): SFB F32‐N18 and BioTechMed-Graz.

References

1) Christoph Aigner, Christian Clason, Armin Rund, Rudolf Stollberger. Efficient high-resolution RF pulse design applied to simultaneous multi-slice excitation. J Magn Reson., 263:33-44, 2016.

2) Armin Rund, Christoph Aigner, Karl Kunisch, Rudolf Stollberger. Simultaneous multislice refocusing via time optimal control. Magn. Reson. Med. 80(4):1416-1428, 2018.

3) Gang Xiao, Renhua Wu, Phillip Zhe Sun. Fast simulation and optimization of pulse-train chemical exchange saturation transfer (CEST) imaging. Phys Med Biol, 60(12):4719-4730, 2015.

4) Jacques Bittoun, Jacques Taquin, Michel D. Sauzade. A computer algorithm for the simulation of any nuclear magnetic resonance (NMR) method. J Magn Reson, 2(2):113-120, 1984.

5) Matt A. Bernstein, Kevin F. King, Xiaohong Joe Zhou. Handbook of MRI Pulse Sequences. Elsevier Academic Press, 2004.


6) Moritz Zaiss and Peter Bachert. Exchange-dependent relaxation in the rotating frame for slow and intermediate exchange - modeling off-resonant spin-lock and chemical exchange saturation transfer. NMR in Biomedicine, 26:5, 2013.

Figures

Figure 1: Evaluation of different Bloch simulations for slice selective inversion. The first column shows the complex RF and slice selective gradient (Gs) as well as the x-, y- and z-components of the magnetization vector calculated with EV and without relaxation. Column 2-4 depict the logarithmic error of the z- and (x,y)-coordinates compared to EV for different Bloch solvers and relaxation times.

Figure 2: Overview of the relative L2 - and L∞ - error compared to the EV results for different Bloch solvers and relaxation times. The last column contains the average runtimes, where the upper five solvers are implemented non-vectorized.

Figure 3: Evaluation of different Bloch-McConnell simulations with the Gaussian pulse. Column 1 shows the RF and the corresponding z-magnetization of the water proton pool calculated with EV. In column 2 and 3, the logarithmic errors of the z- and (x,y)-magnetizations of the water proton pool and the solute proton pool compared to EV are depicted, respectively.

Figure 4: Relative L2 - and L∞ - errors ε of the Bloch-McConnell solvers ASY and SY compared to the EV results for a series of block pulses (BL) and a series of Gaussian pulses (G). The error study is done separately for each coordinate x, y and z and for the water (w) and solute proton pool (s).

Figure 5: Convergence rate of the proposed Bloch-McConnell solvers for a series of RF block pulses (BL) with decreasing time step size τ. The columns show the relative L2 - error in each case compared to EV, and the experimental order of convergence (EOC).

Proc. Intl. Soc. Mag. Reson. Med. 27 (2019)
0760