Investig Magn Reson Imaging. 2019 Mar;23(1):38-45. English.
Published online Mar 28, 2019.
Copyright © 2019 Korean Society of Magnetic Resonance in Medicine (KSMRM)
Original Article

High-Resolution Numerical Simulation of Respiration-Induced Dynamic B0 Shift in the Head in High-Field MRI

So-Hee Lee,1,2 Ji-Seong Barg,1,2 Seok-Jin Yeo,1,2 and Seung-Kyun Lee1,2
    • 1Center for Neuroscience Imaging Research, Institute for Basic Science (IBS), Suwon, Korea.
    • 2Department of Biomedical Engineering, Sungkyunkwan University, Suwon, Korea.
Received August 29, 2018; Revised November 29, 2018; Accepted February 01, 2018.

This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/3.0/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Purpose

To demonstrate the high-resolution numerical simulation of the respiration-induced dynamic B0 shift in the head using generalized susceptibility voxel convolution (gSVC).

Materials and Methods

Previous dynamic B0 simulation research has been limited to low-resolution numerical models due to the large computational demands of conventional Fourier-based B0 calculation methods. Here, we show that a recently-proposed gSVC method can simulate dynamic B0 maps from a realistic breathing human body model with high spatiotemporal resolution in a time-efficient manner. For a human body model, we used the Extended Cardiac And Torso (XCAT) phantom originally developed for computed tomography. The spatial resolution (voxel size) was kept isotropic and varied from 1 to 10 mm. We calculated B0 maps in the brain of the model at 10 equally spaced points in a respiration cycle and analyzed the spatial gradients of each of them. The results were compared with experimental measurements in the literature.

Results

The simulation predicted a maximum temporal variation of the B0 shift in the brain of about 7 Hz at 7T. The magnitudes of the respiration-induced B0 gradient in the x (right/left), y (anterior/posterior), and z (head/feet) directions determined by volumetric linear fitting, were < 0.01 Hz/cm, 0.18 Hz/cm, and 0.26 Hz/cm, respectively. These compared favorably with previous reports. We found that simulation voxel sizes greater than 5 mm can produce unreliable results.

Conclusion

We have presented an efficient simulation framework for respiration-induced B0 variation in the head. The method can be used to predict B0 shifts with high spatiotemporal resolution under different breathing conditions and aid in the design of dynamic B0 compensation strategies.

Keywords
Dynamic B0 shift; Dynamic shim; Respiration; Susceptibility; Brain; 7T

INTRODUCTION

In high-field MRI, tissue susceptibility-induced static field (B0) shift can cause significant imaging artifacts and spectroscopic errors. While numerous methods to mitigate the B0 shift have been proposed (1, 2), dynamically varying B0 shift remains particularly difficult to compensate for. Such changes can be caused by hardware drift (3), subject motion (4), and motion of foreign magnetic objects (5). In particular, the respiration-induced B0 shift in the head has been the focus of several experimental studies due to its importance in functional MRI (6, 7) at ultra-high magnetic fields. Experimentally, the dynamic B0 shift can be measured using gradient-echo phase imaging (8), magnetic field probes (9), or from post-processing of multi-channel images (10). While providing valuable real-life B0 changes, experimental measurements of respiration effects have typically been limited to low spatial and temporal resolutions. Moreover, such studies have usually been conducted with only a small number of volunteers, representing a limited set of anatomies and breathing conditions.

In this respect, accurate numerical simulation of a respiration-induced B0 shift could offer an opportunity to investigate different factors affecting B0, with flexibility and high spatiotemporal resolutions. Marques and Bowtell (11) reported a simulation study using a susceptibility-to-B0 calculation algorithm based on the k-space-discretized Fourier space convolution. While fast and popular, the algorithm is known to require large zero-filled buffers in the computational grid to avoid aliasing artifacts (5, 12). Consequently, the study by Marques and Bowtell (11) used a large (7.6 mm isotropic) voxel, partly undermining the advantages of simulation studies.

In this work, we demonstrate dynamic, 1 mm-isotropic-voxel numerical simulation of respiration-induced B0 shift in the head using a detailed 4D human model. We used a recently published rapid B0 calculation algorithm (generalized susceptibility voxel convolution, gSVC) (12) to achieve high computational speeds while eliminating aliasing artifacts. We obtained B0 shift gradients in the head that were consistent with the experimental data in the literature. We propose that the method presented provides an efficient platform to investigate dynamic B0 shift in high-field neuro MRI.

MATERIALS AND METHODS

Human Body Model

The human body model, incorporating realistic organ motion from breathing, was generated from the Extended Cardiac And Torso (XCAT) model, which was originally used for computed tomography and positron emission tomography for radiation dose and image quality calculations (13, 14, 15). The model has also been used in MRI simulations for its utility in multimodal and interventional imaging research (16, 17). XCAT provides a set of adjustable parameters to define the anatomy and dynamics of the human body, including gender, scaling of the head and lung sizes, segmented tissue properties (in our case, interpreted as magnetic susceptibility χ), turning on/off cardiac and respiratory motions, and depth of breathing. Notably, XCAT lets users change the voxel size to create body models with different voxel resolutions. In our work, we used the following parameters: gender = male, scaling = default, voxel size = 1 to 10 mm (varied, isotropic), susceptibility = 0 (air and lung)/−8.5 (blood)/−11 (bone)/−9 (other tissues) ppm, cardiac motion = off, maximum diaphragm motion = 2.0 cm, and maximum anterior/posterior expansion = 1.2 cm (nominally “normal” breathing). The model parameters were saved in a text file, which was read during the XCAT model generation process executed from within the Matlab (Mathworks, Natick, MA, USA) environment. No user-defined interpolation processes were needed to create models with different voxel sizes.

The respiration cycle was set to be 5 sec, starting from a fully exhaled position, designated as frame 1. A total of 10 frames within the cycle were simulated. The XCAT model is capable of generating voxelated whole-body breathing models at specified time points. From this model, we used the upper torso, from below the diaphragm (lumbar vertebrae 2) to the neck, as the time-varying susceptibility region. The B0 shift was then calculated in the brain mask of the first frame. Table 1 shows the first and last (axial) slice indices corresponding to the grids containing the torso portion and the head portion of the body model, at each voxel size. Figure 1 shows the susceptibility model and the brain mask for two frames (1, fully exhaled and 5, fully-inhaled) at two voxel sizes (1 and 7 mm).

Fig. 1
Sagittal cross sections of the brain mask (dotted line) and the susceptibility model at voxel sizes of 0.1 cm (a) and 0.7 cm (b). Tick mark labels indicate voxels. For each voxel size, frames 1 (left) and 5 (right) are displayed out of the 10 frames in the breathing cycle simulated by the 4D XCAT model. Color bar at the bottom indicates susceptibility in ppm (0 = air and lung, −8.5 = blood, −9 = other tissue, −11 = bone).

Table 1
Susceptibility source (torso) and B0 calculation target (head) positions, computational times, and respiration-induced B0 gradients at different human body model resolutions. The computational times are divided between the dipolar field kernel calculation (t1) and dynamic B0 update (t2)

Once XCAT generated the desired 4D model in a set of binary files, a custom Matlab code read the data for subsequent processing and B0 calculations. For all calculations, we used Matlab version R2017b on a 64 bit Windows PC with Intel Xeon CPU and 64 GB of RAM.

Susceptibility-Induced B0 Calculation

In our problem, the B0 change in the head was caused by susceptibility variation in the upper torso. The problem is a case where the susceptibility source grid (upper torso) is physically separated from the B0 target grid (head). In such a case, use of a conventional k-space-discretized B0 calculation algorithm applied to an all-encompassing grid containing both source and target regions is computationally inefficient (18). In the conventional method, the problem of converting a susceptibility map into a B0 change map is treated as that of image transformation, in which an input “image” (susceptibility) is convoluted with a point dipolar field kernel K(rr) to yield an output (B0 map). This method became popular (11) partly because of the existence of a simple expression for the kernel in the Fourier space,

K(k)13kz2k2[1] 

which makes Fourier-domain convolution calculation easy and elegant. However, Eq. [1] is strictly valid when the spatial-domain kernel K(rr) is defined in an infinite space. When discretizing Eq. [1] in the k-space, as is conventionally done, one implicitly assumes that the convolution is calculated in a finite-sized spatial domain. This results in the susceptibility map being replicated periodically outside the domain, leading to aliasing artifacts in the output B0 map (5). Practically, this is mitigated by adding zero-filled buffers around the computational grid (containing both susceptibility and target B0 regions) to sufficiently push away the susceptibility replicas. Proposals have been made to reduce the necessary computational grid size (and the time) by discretizing the kernel K, not in the k-space but in the real space, and also by separating the susceptibility grid from the target grid (18). However, these methods were limited to cases where the two grids have the same size and the voxels are isotropic, as explained by Lee et al. (12).

In this work, we used the gSVC method (12) to calculate the B0 change in the head from the voxelated susceptibility source in the upper torso at each time point in the breathing cycle. The method requires a modified dipolar field kernel that incorporates an analytical solution of the susceptibility-induced B0 field from a uniformly magnetized voxel (cubical in our case) instead of a point dipole. The kernel in the real space is given by

14πi=12j=12k=121i+j+karctanxiyjzkrijk+c·13[2] 

where xi, yj, zk define the eight corners of the voxel, and rijk=xi2+yj2+zk2 . In the second term, c = 1 if the target point is within the voxel and c = 0, otherwise.

Eq. [2] is more complicated than the point dipole kernel but is readily calculable in Matlab. This approach eliminates aliasing artifacts without the need for grid-increasing buffers. Furthermore, in gSVC, the convolution in the real-space is explicitly formulated for separated input (susceptibility) and output (B0) grids, with generally different sizes. This results in the total computational grid having size N=NS+NT−1 in all 3 Cartesian dimensions, where NS and NT denote the source and target grid sizes, respectively. Once the kernel Eq. [2] is computed for the relevant grid of size N, such a kernel can be re-used in a dynamic B0 calculation, where discrete convolution between the (time-dependent) susceptibility and the (fixed) kernel is done for all time points (12).

In order to examine the effect of coarse-graining the susceptibility model, we computed the dynamic B0 fields at voxel sizes ranging from 1 to 10 mm. For each voxel size, we computed the B0 in the brain (with the same voxel size as the susceptibility model) and compared the results. The computational times for the dipolar field kernel calculation (called t1) and dynamic B0 update (called t2) were recorded at each voxel size.

Analysis

In order to focus our attention on the dynamic effects of breathing, the B0 map of the first frame was subtracted from the subsequent frames. In what follows, we assume that the symbol B0 denotes only this difference, namely the dynamic component of the field offset. The respiration-induced B0 in the head is expected to vary slowly in space, dominated by low-order spatial harmonics. Based on this, and in order to compare our results with previous experimental findings (8, 9, 10), we computed the linear gradients of the dynamic B0 shift maps along three directions: Gx = dB0/dx, Gy = dB0/dy, y is not subscript. (Make it “dy”) and Gz = dB0/dz, z is not subscript. (Make it “dz”) where x, y, z were coordinates in the left-to-right, posterior-to-anterior, and feet-tohead directions, respectively, for a patient in the head-first supine position. The gradients were obtained two different ways. First, the B0 values on each slice normal to a given gradient direction (x, y, z) were averaged. The 2D averaged B0 thus obtained was then plotted against the coordinate in the gradient direction and fitted with a straight line. The slope of the line represented the gradient of B0. This was the method used by Van de Moortele et al. (8). Second, we fitted the simulated 3D B0 maps to a volumetric linear function of the form f(x, y, z)= xGx + yGy + zGz + const to obtain three gradients. The fitting in this case was done by minimizing the mean square of the residuals using the Matlab function ‘fminsearch’.

RESULTS

Magnitude of Temporal Variation

Figure 2 shows the temporal variation of the B0 shift shown on the midsagittal plane at voxel sizes 0.1 cm and 0.7 cm. All shifts are shown relative to the first frame. Therefore, the B0 maps show the change in B0 when inhalation occurs and the lung expands. Figure 2a shows that for the highest-resolution simulation (0.1 cm, presumed most reliable), inhalation generally shifts B0 upward. This can be understood as the result of the lung filling with more air, which is relatively more paramagnetic than the tissue. Figure 2b shows that the low-resolution simulation yields B0 shifts that are substantially different from the high-resolution one. From Figure 1, we believe that this was mainly due to the poor representation of the body organs at this voxel size. For example, partial volume errors could lead to overestimation of the total volume of the tissue in the torso at the inhaled position. This could cause a decrease in the B0 shift in the brain because of the negative susceptibility of the (overestimated) tissue.

Fig. 2
Midsagittal view of the respiration-induced B0 shift in the brain at 7T at 10 time points in a breathing cycle. Time order is from left-to-right. The B0 calculation algorithm gSVC was used. (a) and (b) correspond to voxel sizes of 0.1 cm and 0.7 cm, respectively. Note significant differences between B0 shifts computed at 0.7 cm vs. 0.1 cm resolution.

Figure 3 shows the time dependence of the minimum, maximum, and mean values of the B0 shift in the brain at the 0.1 cm voxel size. Also indicated are one standard deviation of B0 at each time point. Over the entire brain, the mean B0 shift peaked at 2.7 Hz at the fully-inhaled position (frame 5). This is consistent with the experimental data by Van de Moortele et al. (8), who reported an overall brain B0 shift due to respiration of 1.45 to 4 Hz at 7T. The maximum local excursion of B0, in contrast, was 7.2 Hz in our model, which occurred at frame 5 at the bottom of the brain. This value, relative to the main magnetic field of 7T, compares favorably with a recent measurement (10) reporting about 2.8 Hz B0 shift at a similar location at 3T.

Fig. 3
Time dependence of B0 shift statistics in the brain during breathing computed at 0.1 cm voxel size. Maximum B0 shift of 7.2 Hz is observed at 2 s (5th frame).

Spatial Gradient

Spatial gradients of the respiration-induced B0 shift can inform dynamic shimming in high-field brain imaging. Figure 4 shows the slice-averaged B0 values and their fitted line in three directions for the fully-inhaled (frame 5) temporal position. Compared to Table 1 in the study by Van de Moortele et al. (8), our results for Gx, Gz are within the range observed experimentally (Gx = +0.0163 Hz/cm in our case, ranged between −0.035 and +0.025 in (8); Gz = −0.365 Hz/cm in our case, ranged between −0.380 and −0.155 Hz/cm in (8)). In contrast, we found higher simulated respiration-induced Gy (−0.236 Hz/cm) compared to their data (8) (−0.093 to −0.032). This discrepancy can be because, in that study, the measurement was made with a surface coil positioned on the posterior side of the head and may not have covered the frontal brain sufficiently.

Fig. 4
Mean B0 shift gradients in the brain at the 5th frame (fully-inhaled), computed at 0.1 cm voxel resolution, along the x (left-to-right, a), y (posterior-to-anterior, b), z (foot-to-head, c) directions. The blue traces correspond to B0 values averaged over slices normal to each direction. The straight lines indicate linear fitting to the data.

Figure 5 shows the temporal variation of the spatial gradients of the respiration-induced B0, determined by volumetric linear fitting. Note that the gradient values at 2 s (frame 5) are not the same as the slopes of the fitted lines in Figure 4, because of the different fitting methods. The peak gradients (at frame 5) in Figure 5 are +0.008, −0.182, and −0.262 Hz/cm for x, y, z, respectively. The magnitude of these values can be compared to those of Vannesjo et al. (9), who also used volumetric fitting. In Figure 2 of that study, the linear dynamic gradient terms are displayed as the field at 10 cm from the isocenter. Also, the x and y direction definitions were swapped compared to ours. Given these facts, the data in that study are approximately in agreement with our data, showing about 0.2 Hz/cm variation for anterior-posterior and feet-head directions and negligible gradient in the left-right direction. The feet-head gradient (Gz) variation was slightly larger than the anterior-posterior gradient (Gx) in the study by Vannesjo et al. (9), while the two gradient variations had the same sign, which is also consistent with our result.

Fig. 5
Plot of B0 shift gradients as a function of time at 0.1 cm voxel resolution. The gradient values were obtained by volumetric linear fitting at each time point. The maximum gradient amplitudes in the y, z directions are observed at 2 s (5th frame). Gradient in the x-direction shows little change.

Voxel Size Dependence

Table 1 shows that the computational time decreased dramatically as the voxel size increased. Both dipolar field kernel calculation (t1) and dynamic B0 update (t2) times scaled approximately inversely with the voxel volume (plot not shown). Applying gSVC to the physically separated source (torso) and target (head) grids allowed time-efficient computation of the B0 shift at high (0.1 cm) isotropic resolution, taking about 90 seconds (t1 + t2) for the 10 frames.

While reduced spatial resolution greatly speeds up the computation, Figure 6 shows the increased problem (lack of reliability) associated with simulation at large voxel sizes. At voxel sizes larger than 0.3 cm, it was observed that the Gy magnitude surpassed that of Gz, which was not consistent with the literature (8, 9). Furthermore, for voxel sizes greater than 0.5 cm, the B0 gradient values became more erratic, departing significantly from the high-resolution cases.

Fig. 6
Plot of B0 shift gradients at frame 5 relative to frame 1 as a function of the voxel size. The gradient values become increasingly more erratic as the voxel size increases, indicating reduced reliability at large voxel sizes (larger than about 0.5 cm).

DISCUSSION

We have utilized a high-resolution 4D human body model and a computationally efficient B0 calculation method (gSVC) to simulate the respiration-induced B0 shift in the brain. Compared to previous simulation work (11), we implemented up to 7.63 = 439 times smaller voxel volumes, thereby capturing details of organ movements associated with respiration. The human body model used allows for the variation of many anatomical and physiological parameters which can be used for future investigations involving different breathing conditions.

We found that respiration with a moderate breathing depth (diaphragm movement of 2.0 cm) caused an average B0 shift in the head of up to about 7.2 Hz at 7T. The spatial distribution of the induced B0 shift showed mainly a z-directional gradient, accompanied by a significant y-directional gradient in our model. These are in agreement with a previous report (9) based on direct field-probe measurements around the head.

We also examined the accuracy of the results as a function of the model spatial resolution. The results showed that, as the voxel size increased, the spatial gradient of the B0 shift substantially changed for voxels larger than 0.5 cm. This indicates that at such a large voxel size, the breathing anatomy is probably not adequately captured to produce a reliable B0 prediction in the head. While our data are limited to a particular human body model, the observed voxel size dependence provides a guideline for future investigations of respiration effects using different susceptibility models and computational methods.

The present work had several limitations. First, the susceptibility model was simplified. In particular, susceptibility changes in the airway accompanying oxygen concentration change were ignored. We believe that this is justified based on the study (19) which showed that the respiration-induced B0 shift was mainly due to geometry and not air susceptibility changes. Second, we ignored cardiac motion and the associated blood flow. The cardiac motion effect was simulated by Lee et al. (20) and was shown to be much smaller than what we observed here. The effect of blood flow is expected to produce B0 perturbation close to the large vessels and might not significantly affect the global B0 gradient, as was calculated here for the whole brain. Finally, we used only a single specific dynamic human body model (using a specific set of parameters within XCAT). This work was primarily aimed at demonstrating a new numerical simulation platform for calculating the respiration-induced B0 shift, and we defer studies on different body models and parameters to future research. Any findings regarding dependence on the model parameters (such as body mass index) should ultimately be confirmed experimentally.

In conclusion, we showed that the respiration-induced B0 shift in the head can be computed with high spatiotemporal resolution utilizing a detailed dynamic human body model and an efficient B0 calculation algorithm. We have greatly improved the fidelity of the human body model used in the simulation compared to a previous study and confirmed aspects of previous experimental observations. Answering specific questions regarding parametric dependence can be pursued in future research to better guide dynamic shim strategies for ultra-high field neuro MRI.

Acknowledgments

This work was supported by the Institute for Basic Science of the Republic of Korea, under grant IBS-R015-D1. We thank Dr. Yoon-Chul Kim from the Samsung Medical Center for informing us of the XCAT program. We acknowledge support from the KSMRM MR Engineering Research Group.

References

    1. Koch KM, Rothman DL, de Graaf RA. Optimization of static magnetic field homogeneity in the human and animal brain in vivo. Prog Nucl Magn Reson Spectrosc 2009;54:69–96.
    1. Kim PK, Lim JW, Ahn CB. Higher order shimming for ultrafast spiral-scan imaging at 3 tesla MRI system. J Korean Soc Magn Reson Med 2007;11:95–102.
    1. Foerster BU, Tomasi D, Caparelli EC. Magnetic field shift due to mechanical vibration in functional magnetic resonance imaging. Magn Reson Med 2005;54:1261–1267.
    1. Hutton C, Andersson J, Deichmann R, Weiskopf N. Phase informed model for motion and susceptibility. Hum Brain Mapp 2013;34:3086–3100.
    1. Bouwman JG, Bakker CJ. Alias subtraction more efficient than conventional zero-padding in the Fourier-based calculation of the susceptibility induced perturbation of the magnetic field in MR. Magn Reson Med 2012;68:621–630.
    1. Zahneisen B, Asslander J, LeVan P, et al. Quantification and correction of respiration induced dynamic field map changes in fMRI using 3D single shot techniques. Magn Reson Med 2014;71:1093–1102.
    1. Zeller M, Kraus P, Muller A, Bley TA, Kostler H. Respiration impacts phase difference-based field maps in echo planar imaging. Magn Reson Med 2014;72:446–451.
    1. Van de Moortele PF, Pfeuffer J, Glover GH, Ugurbil K, Hu X. Respiration-induced B0 fluctuations and their spatial distribution in the human brain at 7 Tesla. Magn Reson Med 2002;47:888–895.
    1. Vannesjo SJ, Wilm BJ, Duerst Y, et al. Retrospective correction of physiological field fluctuations in high-field brain MRI using concurrent field monitoring. Magn Reson Med 2015;73:1833–1843.
    1. Meineke J, Nielsen T. Data-driven correction of B0-off-resonance fluctuations in gradient-echo MRI; Proceedings of the 26th Annual Meeting of ISMRM; Paris, France. 2018. pp. 1172.
    1. Marques JP, Bowtell R. Application of a Fourier-based method for rapid calculation of field inhomogeneity due to spatial variation of magnetic susceptibility. Concepts Magn Reson Part B 2005;25B:65–78.
    1. Lee SK, Hwang SH, Barg JS, Yeo SJ. Rapid, theoretically artifact-free calculation of static magnetic field induced by voxelated susceptibility distribution in an arbitrary volume of interest. Magn Reson Med 2018;80:2109–2121.
    1. Segars WP, Mahesh M, Beck TJ, Frey EC, Tsui BM. Realistic CT simulation using the 4D XCAT phantom. Med Phys 2008;35:3800–3808.
    1. Silva-Rodriguez J, Tsoumpas C, Dominguez-Prado I, Pardo-Montero J, Ruibal A, Aguiar P. Impact and correction of the bladder uptake on 18 F-FCH PET quantification: a simulation study using the XCAT2 phantom. Phys Med Biol 2016;61:758–773.
    1. Koybasi O, Mishra P, St James S, Lewis JH, Seco J. Simulation of dosimetric consequences of 4D-CT-based motion margin estimation for proton radiotherapy using patient tumor motion data. Phys Med Biol 2014;59:853–867.
    1. Lowther N, Ipsen S, Marsh S, Blanck O, Keall P. Investigation of the XCAT phantom as a validation tool in cardiac MRI tracking algorithms. Phys Med 2018;45:44–51.
    1. Paganelli C, Summers P, Gianoli C, Bellomi M, Baroni G, Riboldi M. A tool for validating MRI-guided strategies: a digital breathing CT/MRI phantom of the abdominal site. Med Biol Eng Comput 2017;55:2001–2014.
    1. Dewal RP, Yang QX. Volume of interest-based fourier transform method for calculation of static magnetic field maps from susceptibility distributions. Magn Reson Med 2016;75:2473–2480.
    1. Raj D, Paley DP, Anderson AW, Kennan RP, Gore JC. A model for susceptibility artefacts from respiration in functional echo-planar magnetic resonance imaging. Phys Med Biol 2000;45:3809–3820.
    1. Lee SK, Barg JS, Yeo SJ. Respiration-induced dynamic B0 shifts in the head: numerical simulation based on generalized susceptibility voxel convolution (gSVC); The 6th International Congress on Magnetic Resonance Imaging (ICMRI); Seoul, Korea. 2018.

Metrics
Share
Figures

1 / 6

Tables

1 / 1

PERMALINK