Three‐dimensional spatially resolved phase graph framework

An open‐source spatially resolved phase graph framework is proposed for simulating arbitrary pulse sequences in the presence of piece‐wise constant gradients with arbitrary orientations in three dimensions. It generalizes the extended phase graph algorithm for analysis of nonperiodic sequences while preserving its efficiency, and is able to estimate the signal modulation in the 3D spatial domain.


| INTRODUCTION
The extended phase graph (EPG) formalism, originally introduced by Hennig 1 and further augmented to support 3d gradients with anisotropic Gaussian diffusion, 2 is a powerful tool for analyzing echoes in the Fourier representation. Gradient-induced dephasing is represented as a shift of the corresponding wave vector. The EPG formalism achieves a high computational performance by adding up magnetization components belonging to the same configuration state (in the EPG terminology) but arising from different coherence pathways. Periodic sequences such as rapid imaging with refocused echoes (eg, turbo spin echo, fast spin echo) 3 or true fast imaging with steady-state precession (ie, balanced SSFP) 4,5 can be conveniently described with discrete evolution periods and state grids representing discrete dephasing increments. These sequences are therefore amenable to analysis with the conventional EPG formalism. However, it remains a challenge to explore the magnetization states at arbitrary time points or calculate echo intensities of aperiodic sequences with a large number of RF pulses or varying gradients. In such cases, the number of configuration states may grow exponentially.
This challenge has partially been addressed by Kiselev 6 for one-dimensional gradient actions by indexing all wave vectors, which are represented as delta functions in the Fourier domain, and merging closely spaced components. The explicit involvement of the exact Green's function of the Bloch-Torrey equation 7 in the presence of a constant gradient made it possible to account for diffusion and global transport. Although the theoretical framework presented by Kiselev 6 could readily be extended to three dimensions, the merging algorithm was only effective along a single spatial direction. Therefore, the efficient implementation of the method remained one-dimensional.
In this work, we have extended this method 6 by (1) implementing support for 3D gradients with arbitrary orientation and adding a novel component-merging algorithm for efficient computations; (2) accounting for anisotropic Gaussian diffusion in three dimensions through the corresponding diffusion tensor; (3) abandoning the treatment of the magnetization evolution in the infinite space and establishing a link between the developed technique to imaging. The latter is achieved by a postsimulation module capable of calculating spatial modulation functions and of pictorial visualization of the magnetization evolution in the 3D spatial domain.
We refer to the proposed framework as 3D spatially resolved phase graph (3D-SPRG). The developed technique with the proposed merging algorithm and 3D visualization was validated for three illustrative example sequences: (1) fast off-resonance calculation for pseudo SSFP (pSSFP) 8 as used in MR fingerprinting, (2) a comparison between quasiisotropic and conventional diffusion-weighting steady-state (DW-SSFP) 9 sequences, and (3) advanced spoiler gradient design for spectroscopic imaging localized with PRESS. 10,11 2 | THEORY

| Evolution of magnetization and simulator
Considering unrestricted anisotropic diffusion, the evolution of magnetization during each interpulse interval is defined by the Bloch-Torrey equation 2,6 as follows: The 3D magnetization vector, (r, t), is the difference of the current magnetization Ψ and that at the thermal equilibrium, with the latter defined as [0, 0, 1] t ; D is the macroscopic effective self-diffusion tensor with the corresponding tensor elements D mn = D nm , and any double-occurring index implies a summation over this index taking values x, y, z ("the Einstein notation"); igrI z describes the dephasing effect of the gradient, in which gr is the dot product of the gradient vector g and position vector r; and I z is the generator of rotations around the z-axis: The terms 1 describe T 2 and T 1 relaxation effects, respectively. Beyond the trivial relaxation terms, Equation (1) has the following propagator (Green's function) 6,7 : where r 0 is the source point, and r is the displacement from the source. This propagator has a form of 3 × 3 matrix acting on the initial 3D magnetization vector. The matrix exponent is calculated, keeping in mind that e i I z is the rotation matrix to the angle α around the z-axis, and I 2 z is diagonal, 6 and the determinant of D equals the product of its eigenvalues.
It is convenient to consider Equation (3) in Fourier space, where the last factor e igr 0 tI z in Equation (3) maps to the wellknown operators that shift the 3D wave vectors of already present magnetization components. The pulse sequence can then be represented as a stack of such operators interlaced with the instantaneous rotation of the entire magnetization by the applied RF pulses. The initial spatially constant magnetization, which is proportional to (2π) 3 (k) in the Fourier domain, develops then to a weighted sum of 3D delta functions: (2) 3 12 I 2 z e igr 0 tI z = r − r 0 , t e igr 0 tI z , where index a labels the ath configuration state, and C a denotes its magnetization vector-valued weighting coefficient. The first exponential factor r − r 0 , t in Equation (3) results in additional multiplication with its Fourier transform, which accounts for the attenuation due to diffusion: For clarity, the omitted subscripts in this equation label the elements of the propagator matrix. Note that this formula has not been explicitly presented before. 6 The remaining relaxation effects are taken into account by multiplying such term with the matrix: According to Figure 1, the procedure outputs pairs of C a and k a after each loop, and finishes when the current RF pulse counter N reaches the total number of RF pulses N rf . Within the entire algorithm, all vectors and coefficients are stored in memory as lists. In addition, no rounding or gridding is applied in any of these analytic calculations.

| Wave number merging
To mitigate the computational burden for nonperiodic pulse sequences, we developed a wave vector-merging algorithm. We merge any two coherence components with close k-vectors, k a and k b , into one new coherence component with the vector where the corresponding However, for a large list of coherence components characterized by their 3D vectors k a , it remains a challenge to find the closest neighbor effectively.
F I G U R E 1 Flowchart of the 3D spatially resolved phase graph simulation framework. Abbreviations: DFT, discrete Fourier transform; Kg, grid merging step To avoid the scaling of computation times as N 2 , where N is the number of wave vectors, we propose a neighbor search algorithm that is based on imposing a rectangular 3D grid with the predefined bin size on the k-vectors. The target bins for k-vectors are calculated on the fly by rounding components of the k-vectors to the integer multiples of the specified merging grid step. The coherence components within the same bin are considered neighbors and are replaced by a single component at the weighted-average k-vector k a ′ , while all k-vectors without neighbors remain unaffected. The proposed wave number-merging algorithm is exemplified in Figure 2 for two frequency dimensions.

| Three-dimensional spatial modulation function
The basic building block for the EPG and the following modifications 1,2,12 is the helix-shaped magnetization modulation induced by applied magnetic field gradients. The final magnetization is the sum of such helix-shaped components with the calculated weights. In a spatially unconstrained region, it is only the component with k = 0 that gives a nonzero contribution to the overall signal, the echo.
We now consider a more realistic case of a finite imaged object defined using its proton density, (r). For clarity, we assume spatially constant D, T 1 , and T 2 . Then, the signal in the image domain takes the following form: when properly normalized, the sum in Equation (7) turns into the inverse Fourier transform of the coefficients C a . We term (r) a spatial modulation function, as it is multiplied with the underlying proton density to form the observed signal. The sum in Equation (7) is to be cropped according to the size of (7) S (r) = (r) ∑ a C a e ik a r = (r) (r) .

F I G U R E 2
The grid merging scheme is illustrated for a 2D k-space with conventional gridding and the proposed method. A, In this example, the field contains five vectors in the list before merging. B, With the defined grid merging step Kg, the field is split into nine bins. In the conventional gridding, all vectors are shifted to the nearest grid point and merged (C), whereas in the proposed method (D) the k-vectors are only altered if the corresponding bin contains more than one component. The result of the grid merging will be, in this case, a list of three k-vectors with corresponding weights: one copied from the original input list and the other two calculated as weighted sums of the respective input entries. The proposed method amounts to fewer and smaller inaccuracies compared with the conventional method, as gridding-induced shifts depend primarily on the distance of two vectors before merging (indicated by boxes in [D]), which is typically much smaller than Kg the actually sampled k-space. Note that the multiplication with the proton density to arrive to signal distribution corresponds to a convolution in k-space; therefore, some extension of the k-space bounding box may need to be considered, depending on the properties (eg, spatial frequency content) of (r). Although strictly speaking this treatment is only valid for samples with homogeneous relaxation and diffusion properties (eg, phantom solutions), it may be extended to multi-compartment objects by considering each compartment separately.

| Software implementation
For efficiency, the core calculation routines of the 3D spatially resolved phase graph method were implemented in C ++ as an extension of the original algorithm proposed by Kiselev,6 whereas for ease of use, the interface for collection of the results and transformation to the image domain were implemented in MATLAB (MathWorks, Natick, MA), as shown in Figure 1.
The present algorithm built on Equation (3) assumes constant gradients between the RF pulses. The tasks for which it is important to account for the time-dependent gradients, it is possible to approximate them in a step-wise manner. This can be achieved trivially by introducing a number of dummy RF pulses with 0° flip angle as separators of multiple constant gradient intervals. The time stepping can be freely chosen, as required by the target application. The grid step for the basic version of the proposed 3D wave number merging procedure is defined by the user, to balance the rounding errors originating from this merging step and the computational burden. Based on the investigation of the effect of the merging grid step on the computation time and simulation errors (presented subsequently), an automatic procedure for finding the optimum gridding step was introduced. It is achieved by running the entire simulation iteratively, reducing the grid step Kg for each iteration until the change between successive simulations is reduced below a user-defined tolerance. In the example in this paper, the tolerance was defined as the RMS error of signal over the object space, but one may redefine it as required for the specific application.
After finishing the magnetization evolution calculation for the entire sequence, the postsimulation module is exploited for (1) collecting all of the results in terms of lists of k-vectors and corresponding coefficients, (2) plotting the signal evolution in 3D k-space, and (3) presenting the signal modulation function in the image domain. To visualize the signal evolution in k-space, the generated k-vectors are color-coded by the corresponding signal magnitude coefficients. The discrete Fourier transform is used to transform the magnetization evolution to the image domain according to Equation (7), which was computationally advantageous compared to gridding followed by a fast Fourier transform, as the number of k-vectors in the sampled k-space region is typically small compared with the number of cells in the 3D k-space.

| Off-resonance simulation
An MR fingerprinting pattern with pSSFP 8 was chosen for frequency-response simulation. This sequence was decoupled to address the sensitivity of SSFP-based MR fingerprinting to magnetic field inhomogeneities, and maintains the spin-echo behavior for variable flip angles and TRs. To fulfill the pSSFP condition, a sequence of TE and TR was calculated based on the flip angle (FA) evolution. To simulate the frequency effect, a static background gradient (G stat ) was presumed, applied along z-direction, giving rise to a range of off-resonant frequencies (2π*[−100 100] rad/s) over a spatial range of 128 mm.
To evaluate the effect of k-space grid merging, the entire sequence was simulated multiple times with the adapted grid step Kg for each iteration. The value of Kg was initialized as 20 rad/m (0.003 mm -1 ), which corresponds to about 1/128 of the k-space extent of a readout for a 1-mm voxel, on the order of the gradient moment needed to advance from one k-space sample to another. For each iteration, Kg was reduced by a factor of 0.2. The iterative simulation stopped when the signal amplitude change between two successive iterations was below 1%. To estimate the residual error of simulation results, a conventional Bloch simulation was exploited additionally as a reference. A virtual white-matter tissue 13 with T 1 = 1084 ms and T 2 = 69 ms was assumed, and a pSSFP sequence with the number of pulses (N rf ) = 100 was simulated on a PC equipped with an Intel i7-4770 3.40-GHz CPU and 16GB of RAM. To test the gridding performance, the number of k-vectors, residual error, and computation time were compared for varying merging grid steps. During the off-resonance response visualization, the frequency resolution of 1 rad/s was used.

| Diffusion-weighted SSFP simulation
A quasi-isotropic DW-SSFP sequence 9 is considered to provide a quasi-isotropic diffusion weighting, such as for a diagnosis of stroke in brain white matter by suppressing anisotropic diffusion effects, and was proposed as an alternative to the conventional DW-SSFP sequence. 14 To investigate the properties of quasi-isotropic DW-SSFP, its frequency response was compared with a conventional DW-SSFP sequence for different diffusion-weighting directions. The sequence diagrams are shown in Figure 4A-C.
Each sequence was simulated twice, and for the second simulation the diffusion-weighting direction was rotated by 45° with respect to the first. Other parameters were FA = 25°, N rf = 100, Kg = 1 rad/m, diffusion gradient amplitude = 23.5 mT/m, diffusion gradient duration = 5 ms, and diffusion µm 2 /ms. 15 As in the previous example, a static background gradient G stat was applied along the z-direction. Diffusivity of 0 along the z-axis was set to eliminate possible diffusion weighting by G stat .

| Point-resolved spectroscopy simulation
The cutting-edge DOTCOPS (dephasing optimization through coherence order pathway selection) method 16 optimizes the crusher gradient scheme of the PRESS sequence to eliminate unwanted coherence pathways. In comparison with the identical spoiler gradients along all three axes used in conventional PRESS, DOTCOPS adopts a different spoiler strategy along the third axis (as shown in Supporting Information Section S1). For illustrative purposes, the signal from the volume-ofinterest edges was assumed where the refocusing FAs deviate from the nominal values (ie, α 1 = α 2 = α 3 = 90°), representing half of the nominal FAs for the refocusing pulses. The following parameters for a PRESS-based spectroscopic imaging sequence were assumed: volume size = 48 mm 3 , matrix = 16 × 16 × 16, TE =30 ms with TE 1 = 14 ms (from the excitation to the first spin echo) and TE 2 = 16 ms (from the first spin echo to the second), k-space merging grid step = 1 rad/m, and spoiling strength for each crusher gradient = 500 m -1 ; relaxation effects were ignored. To simulate the residual frequency dispersion due to shimming 17 imperfections over the entire measurement time, additional static gradients were supposed to be 0.1/−0.2/0.3 Hz/cm on x/y/z axes, respectively. To visualize the magnetization evolution in k-space, four snapshots were taken at time 1 = TE 1 , time 2 = TE, time 3 = TE + 500 ms, and time 4 = TE + 1000 ms, respectively.
The signal acquisition in spectroscopic imaging was assumed to be a purely phase-encoded FID, leading to a cubic-sampled region in k-space. To account for all possible spurious signals received during the entire spectroscopic imaging sequence, a single simulation for the k-space center was performed, and the coherence vectors within the sampled k-space area were included.

| Off-resonance simulation
In the pSSFP frequency-response simulation with 100 RF pulses (Figure 3), the implemented TE, TR, and FAs were derived based on the pSSFP condition ( Figure 3A). When iteratively reducing the merging grid resolution Kg, the number of k-vectors increased ( Figure 3B) for each iteration, but far below the theoretical exponential growth, shown as a dashed black line. The corresponding computation time and residual error are presented in Figure 3C. In contrast to a computation time of about 2.2 hours needed for Kg = 0.032 rad/m, the signal spectra for Kg = 0.8 rad/m and Kg = 0.16 rad/m manifested a comparable RMS error (in comparison to the reference spectra calculated from the conventional Bloch simulation), also in frequency spectra ( Figure 3D). The spectra from the Bloch simulation are not shown.

| Diffusion-weighted SSFP simulation
In the conventional DW-SSFP ( Figure 4D,E, solid lines), the simulation reveals frequency-response profiles, showing similar shapes but different scaling for the two weighting directions (0° and 45°). Within the flat part of the frequency response, the discrepancy between the two weightings remains practically constant within a wide range of off-resonance frequencies, and the signal remains stable along the echo train after the steady state is achieved. In contrast, the spectra calculated for the quasi-isotropic DW-SSFP sequence for the two orientations ( Figure 4D,E, dashed lines) show differences across the entire frequency range. Some similarities in shape of the off-resonance frequency response (beyond ± 85 Hz) can be observed, but the scaling still departs significantly from the expected isotropic behavior. Further discrepancies as well as signal oscillations along the echo train are observed on-resonance (see, for example, echo numbers 49 and 50). These simulations reveal the signal-evolution patterns that depart notably from the perfectly isotropic diffusion weighting assumed in the original publication. 9 Our results suggest that an oscillating steady state of directionally weighted signals is created in quasi-isotropic DW-SSFP, similar to that observed in another study. 18 Figure 5 presents the evolution of magnetization in 3D k-space for the PRESS sequence when considering static gradients (from time 1 to time 4). The complete set of resulting k-vectors is presented in Supporting Information Sections S2 and S3 for a global view, while k-vectors lying inside the cubic sampled area (defined by the matrix/volume size) are presented as a 3D k-space map. In addition, k-vectors inside the sampling area were projected onto the three orthogonal planes and transformed into the image domain. To highlight the relevant k-vectors, the vectors with weighting coefficients C below a certain user-defined threshold (10 -8 ) were not plotted.

GAO et Al.
When the conventional spoiler scheme is used ( Figure 5A), not only the desired k-vector but also spurious k-components are generated at the TE (t 2 ) and produce artifacts on the three projection planes for both magnitude ( Figure 5B) and phase ( Figure 5C) images after the discrete Fourier transform. These artifacts manifest themselves as a spatial modulation of the signal visualized for each plane, based on the projection of the corresponding k-vectors ( Figure 5A at t 2 ).
Due to the static gradients assumed during the acquisition, the k-vectors evolve (denoted by "*" as shifting from the k-space center) between time 3 and time 4, but the corresponding weighting coefficients stay constant, which produces additional phase-only artifacts ( Figure 5C at t 3 and t 4 ). These phase artifacts are different for the three projection planes due to the distinct static gradients for the three orthogonal axes. Consequently, this is expected to give rise to artifacts in the acquired spectra if reconstructed without proper phase correction. Note that the spurious k-vector at TE only arises from the spatial-selection gradients, as the static gradients are fully balanced at this point.
When the DOTCOPS scheme is used, no echo is observed at t 1 . The spurious signals are also canceled at TE (t 2 ), and only the double spin echo (k = [0, 0, 0] t ) remains ( Figure 5D at t 2 ). Although the magnitude images are not affected at t 3 and t 4 ( Figure 5E), phase artifacts are generated ( Figure 5F) when additional static gradients are considered. Such phase artifacts may propagate into spectroscopic analysis and quantification if no appropriate phase correction is applied. This F I G U R E 3 The frequency-response simulation of the pseudo-SSFP (pSSFP) sequence is presented for iterative adaptation of Kg. The result from a conventional Bloch simulation was set as a reference. For each iteration, a sequence simulation with 100 RF pulses was performed. A, The implemented TEs, TRs, and flip angles (FAs) were calculated based on the pSSFP condition in the protocol parameters graph. B, In the burden graph, the number of k-vectors is displayed as a function of the number of pulses, where the dashed black line depicts the exponentially growing number of k-vectors without the merging algorithm. For each iteration, the computation time and the residual errors compared with the reference are presented for different Kg in the efficiency graph (C), and the frequency-response spectra after 100 pulses are depicted in (D) underlines the importance of careful shimming in spectroscopic imaging, as shim imperfections may cause signal modulation during volume selection.

| DISCUSSION
In this paper, we propose an integrated simulation framework, 3D spatially resolved phase graph, for simulating echo amplitudes in arbitrary pulse sequences efficiently and intuitively. It extends the core algorithm of Kiselev 6 by adding support for arbitrary gradients and anisotropic diffusion and implementing a new k-merging algorithm for an efficient computation. Additionally, a MATLAB-based visualization module is supplied for pictorial visualization of spin evolutions in the 3D spatial domain.
The present method is the solution to the Bloch-Torrey equation using the algorithm of characteristics 2 without relying on any periodicity requirements in the pulse sequence, in contrast to the classical EPG. 12 Meanwhile, our method preserves other advantages of the EPG, such as providing an analytical solution for the time intervals of constant gradients and maintaining a finite number of the dephasing components. Compared with the requirement of the infinitesimal time interval and infinitesimal slicing of space for solving the Bloch-Torrey equation directly, this leads to favorable computation times.
As for other algorithms, 19 there is a restriction of using hard RF pulses. This can be overcome by splitting the long pulse into multiple low-FA pulses. In this way, when necessary, the present method can be augmented by considering involving infinitesimal time steps for some parts of the sequence. The same trick can be applied to time-varying gradients as well, which is represented in step-wise periods of constant gradients, separated by 0° RF pulses. Such calculation-demanding sequences can be handled effectively by the proposed algorithm and the high performance of the numerical implementation.
Three illustrative examples are presented to demonstrate visualization options of the developed tool for better understanding of the signal behavior. The grid-merging algorithm used is validated for reducing calculation time efficiently, while preserving the desired accuracy. A pictorial 3D presentation leads to a better understanding the underlying spin physics, which cannot be visualized with the one-dimensional algorithms 6,12 as proposed previously. In addition to visualization, because the postsimulation step is decoupled from the simulation kernel, the entire off-resonance signal spectra with an arbitrary frequency resolution can be recalculated based on a single magnetization simulation without repeating the entire simulation procedures.
The proposed simulation algorithm provides a pictorial representation of the signal evolution and a flexible 3D sequence simulation framework. Beyond our simple examples, the proposed approach will be beneficial in more advanced applications such as compressed-sensing trajectory optimization, by characterizing eddy current-induced time-varying background gradients, or intra-echo-train prospective F I G U R E 4 The pulse-sequence diagram of SSFP with background gradient is presented in (A), where the diffusion effect is encoded by conventional (B) or quasi-isotropic (C) diffusion-weighting gradients. Each sequence is calculated twice, with the diffusion-weighting gradients rotated in place by 45° in the second simulation compared with the first. D,E, The corresponding frequency spectra after 49 pulses and 50 pulses are depicted, respectively, in which c 0 and c 45 refer to the conventional diffusion-weighted SSFP simulations, and qi 0 and qi 45 refer to the quasiisotropic simulations