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¹, 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² and GPU acceleration³. 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₀ 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⁴ 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.