Fast simulation of off-resonance artifacts in MRI using FORECAST (Fourier-based Off-REsonanCe Artifact Simulation in the STeady-State)

Frank Zijlstra^{1}, Job G Bouwman^{1}, Ieva Braškutė^{1}, and Peter R Seevinck^{1}

Numerical
simulation of MRI is an important investigative tool for MR research. Examples
of its application include pulse sequence design, simulation of various
artifacts, such as distortions around metal^{1}, and testing of
reconstruction algorithms. Most MRI simulators are based on the discrete time
solution of the Bloch equations, commonly referred to as Bloch simulations. As
simulation of a large number of spins is computationally intensive, Bloch
simulations have been accelerated with parallel computing^{2} and GPU
acceleration^{3}. However, the simulation size still scales quadratically (2D)
or cubically (3D) with increasing field of view or resolution.

In this abstract we demonstrate an efficient method to accurately simulate off-resonance artifacts with a lower computational complexity than Bloch simulations, as long as some assumptions about the pulse sequence can be made.

*Theory:*

In Cartesian 2D
gradient echo imaging with perfect encoding gradients, the basic signal
equation in the presence of B_{0} inhomogeneity is given by:

$$s(k_x,k_y) = \int M_\perp(x,y,z,t',t) \cdot e^{-i2\pi(k_x x+k_y y+\gamma\kern-0.4em- \Delta B_0(x,y,z)t'))} dx dy dz$$

Here, $$$M_\perp$$$ is the transverse signal evolution without encoding, $$$t$$$ is the time since the start of the pulse sequence, and $$$t’$$$ is the time since the last excitation (i.e. for the $$$n$$$'th repetition: $$$t' = t - (n-1) {TR}$$$). Note that $$$k_x = \gamma\kern-0.8em- (t' - TE) G_x$$$.

Assume that $$$M_\perp$$$ does not depend on $$$t$$$, i.e. the signal evolution without encoding is the same during every repetition, which primarily requires a repeating pulse sequence with the MR signal in a steady-state. In this case the signal equation is almost exactly the Fourier transform of $$$M_\perp$$$, with the exception of effects depending on $$$t'$$$ (e.g. $$$\Delta B_0$$$). In this situation, the Bloch simulation is closely related to a naïve evaluation of the 2D Discrete Fourier Transform.

We propose to evaluate this signal equation using Fast Fourier Transforms of $$$M_\perp(t')$$$ for each $$$t'$$$ during data sampling (i.e. one FFT per readout sample). Factors that do not depend on $$$t$$$ can be encoded in $$$M_\perp(t')$$$ , for example: selective excitation, T1 weighting, and T2 decay.

*Experimental validation:*

We acquired a fast 3D gradient spoiled SSFP scan of a titanium cylinder placed in agar gel (matrix 256x256x90, resolution 1 mm, TE/TR 2.1/7.1, flip angle 30°) at 1.5T (Achieva, Philips, Best, The Netherlands).

We simulated a 2D slice through the center of
the cylinder with both Bloch simulations and FORECAST. An analytical model of
the cylinder was used to create a spin density and
$$$\Delta B_0$$$-map^{4}
with 512x512
spins (i.e. 4 spins per voxel). For the Bloch simulation the actual RF pulses
were simulated, while FORECAST only used the RF bandwidth to suppress
signal from spins that were not excited. For simplicity, no slice excitation
was simulated. We also simulated the full 3D scan with 512x512x180 spins (i.e.
8 spins per voxel) using FORECAST.

Finally we compared the computation times of both methods by performing 2D simulations with matrix sizes varying from 4x4 to 256x256.

Figure 1 shows the acquired scan and 2D simulated images. FORECAST shows only minor differences compared with Bloch simulation, and both methods correspond well to the actual scan except for the field induced by the cuboid-shaped agar phantom, which was not simulated. Simulation time was 2478 seconds for Bloch simulations and 6.5 seconds for our method, on a single core of a Xeon E5-1607 CPU. For the 3D simulation, the simulation time using FORECAST was 883 seconds.

Figure 2 shows the computation times of 2D simulations with increasing matrix size; the acceleration is linear with the matrix size.

We demonstrated FORECAST, a fast alternative to Bloch simulations for the simulation of artifacts resulting from off-resonance in steady-state pulse sequences. For 2D simulation, acceleration factors up to over 350 times were achieved. Numerical results confirm that the acceleration factor increases linearly with the simulation size. Theoretically, for 3D simulation we expect a quadratic acceleration.

In theory, our method generalizes to various pulse sequences, as long as the aforementioned requirements are met. For example, simulation of spin echo and non-Cartesian trajectories is possible. A hybrid approach including traditional Bloch simulation could be used to simulate additional effects, such as the transient state at the start of an acquisition, or RF excitation profiles. Further evaluation with different types of pulse sequences is warranted to test the method and explore its limits.

We foresee a range of applications that may benefit from FORECAST, including estimation of encoding artifacts around metal, QSM, and background field removal algorithms. The unprecedented speed of FORECAST may pave the way for new applications.

1. Bakker CJG, Bhagwandien R, Moerland MA, Ramos LMP. Simulation of susceptibility artifacts in 2D and 3D fourier transform spin-echo and gradient-echo magnetic resonance imaging. Magnetic Resonance Imaging 1994;12:767–774.

2. Stöcker T, Vahedipour K, Pflugfelder D, Shah NJ. High-performance computing MRI simulations. Magn. Reson. Med. 2010;64:186–193.

3. Liu F, Kijowski R, Block WF. MRiLab: Performing Fast 3D Parallel MRI Numerical Simulation on A Simple PC. ISMRM 21st Annual Meeting 2013, Salt Lake City, Utah, USA, April 20-26, 2013.

4. Bouwman JG, Bakker CJG. 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.Figure 1. Cylinder
phantom SSFP scan and corresponding simulation using Bloch simulation and FORECAST.

Figure 2. Graphs showing 2D simulation time (A) of
Bloch simulation and FORECAST and the corresponding acceleration factor from
Bloch simulation to FORECAST (B) for increasing simulation size (from 4x4 to
256x256).

Proc. Intl. Soc. Mag. Reson. Med. 24 (2016)

1936