Simulation of rapid MR sequences is often necessary for characterising transient and steady-state signal behaviour, a tool increasingly used for obtaining accuracy in quantitative imaging. There are broadly two approaches: i) isochromat-summation and ii) Extended Phase Graph¹ (EPG) simulations. Of these, isochromat-summation is more intuitive but its accuracy depends on the number and positions of simulated isochromats. By exploring the link with EPG we derive criteria for ensuring isochromat-based simulations are accurate. Matlab code to reproduce all results may be downloaded from https://github.com/mriphysics/EPG-isochromat.

Consider a voxel with uniform magnetization density and relaxation properties. Rotation due to RF pulses is constant within this voxel and is characterized by flip angle θ and phase angle φ. In order to simulate generation of echoes the voxel is further subdivided into a distribution of ‘isochromats’ at positions r each characterised by gradient dephasing angle ψ, where

$$\psi=\gamma\int_0^\tau\mathbf{G}(t')\cdot\mathbf{r}dt'=\mathbf{\Delta\,k}\cdot\mathbf{r}.\qquad\text{[1]}$$

Bloch simulations are performed, and the predicted signal $$$\tilde{S}$$$ is computed by averaging transverse magnetisation M₊ over all isochromats:

$$M_+(\psi)=M_x(\psi)+iM_y(\psi)$$

$$\tilde{S}=\int_0^{2\pi}M_+(\psi^\prime)d\psi^\prime\approx\frac{1}{N}\sum_{m=0}^{N-1}M_+(\psi_m)\qquad\text{[2]}$$

The integral is approximated as a sum over a set of angles $$$\psi_m$$$. Choice of sampling scheme is the key issue for isochromat-summation simulations. The EPG method approaches the same problem in the reciprocal (Fourier) domain, by representing magnetization in terms of configuration states:

$$\tilde{F}_+(\mathbf{k})=\int_V{M_+\,e^{-i\mathbf{k}\cdot\mathbf{r}}}d^3\mathbf{r}\qquad$$

$$\tilde{F}_-(\mathbf{k})=\int_V{M_-\,e^{-i\mathbf{k}\cdot\mathbf{r}}}d^3\mathbf{r}\qquad$$

$$\tilde{Z}(\mathbf{k})=\int_V{M_z\,e^{-i\mathbf{k}\cdot\mathbf{r}}}d^3\mathbf{r}\qquad\text{[3]}$$

For a sequence with regular timing (fixed Δk), the possible k values generated form a set of discrete values:

$$\mathbf{k}=n\mathbf{\Delta k}\qquad n\in\mathbb{Z}\qquad{[4]}$$

such that the magnetization distribution can be written as a sum over integer-indexed states:

$$M_+(\psi)=\sum_{n=-\infty}^{\infty}{\tilde{F}_n e^{in\psi}}\qquad$$

$$M_-(\psi)=\sum_{n=-\infty}^{\infty}{\tilde{F}^*_{-n} e^{in\psi}}\qquad$$

$$M_z(\psi)=Re\left\lbrace\sum_{n=0}^{\infty}{\tilde{Z}_n e^{in\psi}}\right\rbrace\qquad\text{[5]}$$

Numerical predictions are obtained by following a series of algorithmic rules^2. Once numerical values for the coefficients $$$\tilde{F_n}$$$, $$$\tilde{F_{-n}}$$$ and $$$\tilde{Z_n}$$$ are known, we may determine M₊(ψ) for any ψ. A key general assumption is that states with k≠0 are fully dephased: only the k=0 state contributes measurable signal. By comparison with Eq.2 we identify $$$\tilde{S}=\tilde{F}_0$$$.

The DFT may be represented as follows:

$$\tilde{X}_n=\sum_{m=0}^{N-1}{x_m e^{-2\pi\,i\,n\frac{m}{N}}}\quad\,n,m=0,\dots,N-1\qquad\text{[6]}$$

where N is the number of isochromats in the ensemble. If the isochromat angles are selected according to

$$\psi_m=\frac{2\pi\,m}{N}\quad\,m=0,\dots,N-1\qquad\text{[7]}$$

then:

$$\text{DFT}(M_+)=\sum_{m=0}^{N-1}{M_+(\psi_m)e^{-2\pi\,i\,n\frac{m}{N}}}$$

$$\qquad=\sum_{m=0}^{N-1}{M_+(\psi_m) e^{-i\,n\psi_m}}$$

$$\qquad\qquad\qquad=\tilde{F}_n\quad\quad\,n=-\left\lfloor\frac{N}{2}\right\rfloor,\dots,\left\lfloor\frac{N-1}{2}\right\rfloor\text{[8]}$$

Hence, the DFT may be used to recover the EPG representation directly from an isochromat simulation, with the resulting Fourier coefficients corresponding to the configuration state coefficients. There are two restrictions: i)$$$\psi_m$$$ must be selected according to Eq.7; ii) N≥K where K is the number of non-zero transverse configuration states– if N<K then aliasing will occur (i.e. Nyquist sampling condition).

Maximum configuration order of selected sequences

Figure 1 depicts phase diagrams for SPGR (Fig.1a) and FSE (Fig.1b) sequences. For SPGR we are interested in F₀ directly after each RF pulse; after p pulses K_SPGR=2p-1. For FSE sequences we are interested in F₀ mid-way between refocusing pulses; after p refocusing pulses K_FSE=4p-1.

1. SPGR with differing numbers of Isochromats

SPGR was simulated for p=100 RF pulses. Figure 2 demonstrates exact agreement between EPG and isochromat simulations for N=199(=K_SPGR). Figure 3 shows the same results computed with smaller N: for N<199, FFT(M+) shows aliasing effects. These manifest themselves as errors in $$$\tilde{S}$$$ for N<100, when aliasing corrupts the prediction at n=0.

2. FSE with differing numbers of Isochromats

FSE was simulated for p=100 refocusing pulses. Figure 4 shows that exact agreement is obtained for N=399(=K_FSE). Aliasing results when N<K_FSE however errors in $$$\tilde{S}$$$ only result when N≤200 and aliases corrupt the n=0 prediction. There is a marked difference in error for odd and even N.

3. SPGR with diffusion

SPGR was simulated for p=1500 pulses with N=100 for a range of Φ, with and without diffusion implemented by circular convolution³ (Figure 5). Without diffusion, N=100 results in large errors in the isochromat calculation; inclusion of diffusion drastically reduces errors since it acts as a low-pass filter on F_n – i.e. reduced numbers of isochromats can still reproduce the isochromat signal distributions without aliasing.

We demonstrate that isochromat-simulations are equivalent to EPG by simple use of DFT when dephasing angles ψ are appropriately selected. As a result EPG predictions may be deduced from isochromat-summations and vice versa. EPG indicates that isochromat-summation with N isochromats and p pulses can only guarantee accuracy if N≥p for SPGR and N≥2p for FSE. If N is reduced, errors occur in a predictable way due to aliasing. However diffusion effects can be used to maintain accuracy with smaller numbers of isochromats (N<<p) by acting as low pass filter in the Fourier (EPG) domain. Insights gained may pave the way for faster simulations by using a mixture of domains to increase computational efficiency.