Extended phase graph

In magnetic resonance imaging and nuclear magnetic resonance the extended phase graph (EPG) is a mathematical framework used to track how magnetization evolves in a voxel through a series of radiofrequency (RF) pulses and gradients. While traditional Bloch simulations track individual spins in the spatial domain, EPG operates in the Fourier domain.^[2] The magnetization vectors along one dimension of a voxel are stored as Fourier coefficients. Each Fourier coefficient represents the weighting of a helix that twists along one dimension of the voxel with an integer number of twists. Performing a weighted sum of all the helices gives the magnetization vectors in the spatial domain. In this representation, a gradient along one dimension of the voxel that causes an integer number of twists simply shifts the index of each Fourier coefficient by that integer, making EPG very computationally efficient at modeling sequences with many RF pulses and gradients. EPG was first proposed by Jürgen Hennig to model multiecho sequences^[2] and is now used to model magnetic resonance fingerprinting,^[3] turbo spin echos,^[4] and other sequences.^[5]

Definition

In the Bloch equations, within a voxel each spin along the dimension $z$ can be described by a magnetization vector $[M_{x}(z),M_{y}(z),M_{z}(z)]^{T}$ where every element is real valued.^[6] But $M_{x}(z)$ and $M_{y}(z)$ can also be described by the complex number $M_{+}(z)=M_{x}(z)+iM_{y}(z)$ .^[7] Using a similarity matrix $S$ we can convert between the real valued basis and the complex transverse magnetization basis:

${\begin{aligned}&S={\begin{bmatrix}1&i&0\\1&-i&0\\0&0&1\end{bmatrix}},&\mathbf {m} (z)={\begin{bmatrix}M_{+}(z)\\M_{-}(z)\\M_{z}(z)\end{bmatrix}}=S{\begin{bmatrix}M_{x}(z)\\M_{y}(z)\\M_{z}(z)\end{bmatrix}}\\\end{aligned}}$

where $M_{-}$ is the complex conjugate of $M_{+}$ . To perform a rotation on the complex magnetization vector $\mathbf {m}$ due to an RF pulse with angle $\alpha$ and phase $\phi$ , we first apply a change of basis with $S^{-1}$ , then apply standard cartesian rotation matrices $R_{z}(\phi )$ , $R_{x}(\alpha )$ ,^[8] then another change of basis with $S$ :^[9]

${\begin{aligned}&R_{z}(\phi )={\begin{bmatrix}\cos \phi &-\sin \phi &0\\\sin \phi &\cos \phi &0\\0&0&1\end{bmatrix}},\;\;R_{x}(\alpha )={\begin{bmatrix}1&0&0\\0&\cos \alpha &-\sin \alpha \\0&\sin \alpha &\cos \alpha \end{bmatrix}}\\&\\&T_{\phi }(\alpha )=SR_{z}(\phi )R_{x}(\alpha )R_{z}(-\phi )S^{-1}\\&\\&T_{\phi }(\alpha )={\begin{bmatrix}\cos ^{2}{\frac {\alpha }{2}}&e^{2i\phi }\sin ^{2}{\frac {\alpha }{2}}&-e^{i\phi }\sin \alpha \\e^{-2i\phi }\sin ^{2}{\frac {\alpha }{2}}&\cos ^{2}{\frac {\alpha }{2}}&ie^{-i\phi }\sin \alpha \\-{\frac {i}{2}}e^{-i\phi }\sin \alpha &{\frac {i}{2}}e^{i\phi }\sin \alpha &\cos \alpha \end{bmatrix}}\\&\\&\mathbf {m} ^{+}(z)=T_{\phi }(\alpha )\mathbf {m} (z)\end{aligned}}$

Where $\mathbf {m} ^{+}$ denotes the magnetization vector after the rotation. In the EPG representation we perform a Fourier decomposition of each element of the complex magnetization vector:^[2]

${\begin{aligned}&M_{+}(z)=\sum _{k=-\infty }^{\infty }F_{+}(k)e^{i2\pi kz}\;\Leftrightarrow \;F_{+}(k)=\int _{-0.5}^{0.5}M_{+}(z)e^{-i2\pi kz}dz\\&M_{-}(z)=\sum _{k=-\infty }^{\infty }F_{-}(k)e^{i2\pi kz}\;\Leftrightarrow \;F_{-}(k)=\int _{-0.5}^{0.5}M_{-}(z)e^{-i2\pi kz}dz\\&M_{z}(z)=\sum _{k=-\infty }^{\infty }Z(k)e^{i2\pi kz}\;\Leftrightarrow \;Z(k)=\int _{-0.5}^{0.5}M_{z}(z)e^{-i2\pi kz}dz\\\end{aligned}}$

If $F_{+}(k)=0$ , $F_{-}(k)=0$ , and $Z(k)=0$ when $k>n$ or $k<-n$ . We can write the summations as being from $k=-n$ to $n$ . We can also write the Fourier coefficients as a vector valued function:

${\begin{aligned}&\mathbf {f} (k)={\begin{bmatrix}F_{+}(k)\\F_{-}(k)\\Z(k)\end{bmatrix}},&\mathbf {m} (z)=\sum _{k=-n}^{n}\mathbf {f} (k)e^{i2\pi kz}\\\end{aligned}}$

RF pulses

When performing a rotation, because $T_{\phi }(\alpha )$ is constant across $k$ it can enter the sum:

${\begin{aligned}&\mathbf {m} ^{+}(z)=T_{\phi }(\alpha )\mathbf {m} (z)=\sum _{k}T_{\phi }(\alpha )\mathbf {f} (k)e^{i2\pi kz}\\&\mathbf {f} ^{+}(k)=T_{\phi }(\alpha )\mathbf {f} (k)\end{aligned}}$

It can be seen in the equation above that the rotation matrix can be applied directly to the Fourier coefficients to produce updated Fourier coefficients.^[1]

Gradients

A gradient in the z dimension that creates one twist from $z=-0.5$ to $z=0.5$ can be applied with the matrix:

${\begin{aligned}&T(z)=SR_{z}(2\pi z)S^{-1}={\begin{bmatrix}e^{i2\pi z}&0&0\\0&e^{-i2\pi z}&0\\0&0&1\end{bmatrix}}\\\end{aligned}}$

This matrix is independent of $k$ so it can enter the Fourier sum:^[10]

${\begin{aligned}&\mathbf {m} ^{+}(z)=T(z)\mathbf {m} (z)=T(z)\sum _{k}\mathbf {f} (k)e^{i2\pi kz}=\sum _{k}T(z)\mathbf {f} (k)e^{i2\pi kz}\\\end{aligned}}$

Looking at each element of $\mathbf {m} ^{+}$ :

${\begin{aligned}&M_{+}^{+}(z)=\sum _{k}e^{i2\pi z}F_{+}(k)e^{i2\pi kz}=\sum _{k}F_{+}(k)e^{i2\pi (k+1)z}\\&M_{-}^{+}(z)=\sum _{k}e^{-i2\pi z}F_{-}(k)e^{i2\pi kz}=\sum _{k}F_{-}(k)e^{i2\pi (k-1)z}\\&M_{z}^{+}(z)=M_{z}(z)\\\end{aligned}}$

Performing a substitution with ${\hat {k}}=k+1$ for the first element we get:

${\begin{aligned}&M_{+}^{+}(z)=\sum _{{\hat {k}}=-n+1}^{n+1}F_{+}({\hat {k}}-1)e^{i2\pi {\hat {k}}z}\\\end{aligned}}$

If $F_{+}(n)=0$ and $F_{+}(-n-1)=0$ the above summation is equal to:^[11]

${\begin{aligned}&M_{+}^{+}(z)=\sum _{{\hat {k}}=-n}^{n}F_{+}({\hat {k}}-1)e^{i2\pi {\hat {k}}z}\\\end{aligned}}$

and it can be seen that $F_{+}^{+}(k)=F_{+}(k-1)$ . By the same logic $F_{-}^{+}(k)=F_{-}(k+1)$ . So gradients in the EPG formalism simply shift the indices of the Fourier coefficients.^[12]^[1]

Relaxation

The relaxation and recovery matrices given by the Bloch equations for a real valued magnetization vector $\mathbf {v}$ are:^[7] ${\begin{aligned}&E={\begin{bmatrix}e^{-t/T2}&0&0\\0&e^{-t/T2}&0\\0&0&e^{-t/T1}\end{bmatrix}},&\mathbf {r} ={\begin{bmatrix}0\\0\\M_{0}(1-e^{-t/T1})\end{bmatrix}}\\&\\&\mathbf {v} ^{+}=E\mathbf {v} +\mathbf {r} \end{aligned}}$

Substituting with $\mathbf {v} =S^{-1}\mathbf {m}$ :

${\begin{aligned}&S^{-1}\mathbf {m} ^{+}=ES^{-1}\mathbf {m} +\mathbf {r} \\&\mathbf {m} ^{+}=SES^{-1}\mathbf {m} +S\mathbf {r} \end{aligned}}$

$S$ and $E$ are commutative so $SES^{-1}=ESS^{-1}=E$ . Also $S\mathbf {r} =\mathbf {r}$ , so:

${\begin{aligned}&\mathbf {m} ^{+}(z)=E\mathbf {m} (z)+\mathbf {r} =E\sum _{k}\mathbf {f} (k)e^{i2\pi kz}+\mathbf {r} =\sum _{k}E\mathbf {f} (k)e^{i2\pi kz}+\mathbf {r} \\\end{aligned}}$

$\mathbf {r}$ is a constant vector but we can still Fourier decompose it to $\mathbf {r} (z)=\sum _{k}\mathbf {g} (k)e^{i2\pi kz}$ where $\mathbf {g} (k)=\mathbf {r}$ if $\mathbf {k} =0$ else $\mathbf {g} (k)=\mathbf {0}$ , so:^[13]

${\begin{aligned}&\mathbf {m} ^{+}(z)=\sum _{k}E\mathbf {f} (k)e^{i2\pi kz}+\sum _{k}\mathbf {g} (k)e^{i2\pi kz}=\sum _{k}(E\mathbf {f} (k)+\mathbf {g} (k))e^{i2\pi kz}\\\end{aligned}}$

So it can be seen that to apply a relaxation, you can multiply every Fourier coefficient with the relaxation matrix $E$ and then add the recovery vector $\mathbf {r}$ to the 0th coefficient.^[14]

Total signal in voxel

The total transverse signal in the voxel is obtained by integrating over $M_{+}(z)$ :

${\begin{aligned}&M_{+total}=\int _{-0.5}^{0.5}\sum _{k=-n}^{n}F_{+}(k)e^{i2\pi kz}\;dz=\sum _{k=-n}^{n}(F_{+}(k)\int _{-0.5}^{0.5}e^{i2\pi kz}dz)=\sum _{k=-n}^{n}F_{+}(k){\frac {1}{i2\pi k}}(e^{i\pi k}-e^{-i\pi k})\\\end{aligned}}$

so when $k$ is nonzero, $F_{+}(k)$ is multiplied by zero and doesn't end up in the sum. If $k=0$ the inner integral is 1. So the total transverse signal in the voxel is just $F_{+}(0)$ . By the same logic the total longitudinal signal is $Z(0)$ .^[11]

Coefficient redundancy

It can be proven that $F_{+}(k)=F_{-}(-k)^{*}$ :^[13]^[15]

${\begin{aligned}&F_{-}(k)=\int M_{+}(z)^{*}e^{-i2\pi kz}dz\\&F_{-}(-k)=\int M_{+}(z)^{*}e^{i2\pi kz}dz=\int (M_{+}(z)e^{-i2\pi kz})^{*}dz=F_{+}(k)^{*}\\\end{aligned}}$

and so in practice only coefficients with $k\geq 0$ are stored and if a gradient occurs that would shift a coefficient to having a negative index, instead that coefficient is conjugated and stored at a non-negative index.

Computation

In practice a discrete number of Fourier coefficients are tracked. A matrix can be set up with $n+1$ frequencies tracked like so:^[14]

${\begin{bmatrix}F_{+}(0)&F_{+}(1)&...&F_{+}(n)\\F_{-}(0)&F_{-}(1)&...&F_{-}(n)\\Z(0)&Z(1)&...&Z(n)\\\end{bmatrix}}$

This matrix can be multiplied by the RF pulse rotation matrix $T_{\phi }(\alpha )$ , or multiplied by the relaxation matrix $E$ and then have ${\textbf {r}}$ added to the first column, to update the coefficients. A gradient that adds one twist along the z dimension will result in the first row being shifted right by one, and the second row being shifted left. Then $F_{+}(0)$ can be set to $F_{-}(0)^{*}$ .^[16]