Quantitative susceptibility mapping (QSM) is a novel imaging technique in magnetic resonance imaging (MRI) developed near 2008. It enables the depiction of tissue’s magnetic property and has found many clinical applications such as detecting micro-bleed, lesions, calcification, etc. In principle, it has potential utility in any diseases with metabolic abnormality in magnetic elements, e.g., iron, calcium, magnesium.

In this post, I will summarize the math behind it.

signal equation with Lorentz sphere correction

In the case of magnetostatics without current source, the Gauss’s Law and Ampere’s Law are given by

\[\begin{align} \nabla \cdot \mathbf{B} =& 0 \\ \nabla \times \mathbf{H} =& 0. \end{align}\]

Note here both \(\mathbf H\) and \(\mathbf B\) are macroscopic and they are related by \begin{align} \mathbf{H} = \mathbf{B}/\mu_0-\mathbf{M} \end{align} where \(\mu_0\) is the magnetic permeability of vacuum and \(\mathbf M\) is the macroscopic magnetization or magnetic moment density. In this note we use bold font to represent vector and unbold font to represent the z-component only.

For isotropic diamagnetic and paramagnetic substances, a linear constitutive relationship between \(\mathbf H\) and \(\mathbf B\) holds \begin{align} \mathbf{B} = \mu\mathbf{H} \end{align} where \(\mu\) is the substance permeability. By writing \(\mu=(1+\chi)\mu_0\), where \(\chi\) is the magnetic susceptibility, we get \(\mathbf{M}=\chi\mathbf{H}\).

With the Lorentz sphere correction, the measured NMR signal is given by \begin{align} \mathbf{B_\text{NMR}}=\mathbf{B}-\frac{2\mu_0}{3}\mathbf{M} =\left(1+\frac{\chi}{3}\right)\mu_0 \mathbf{H}. \end{align}

Taking the curl of Ampere’s law, we get \begin{align} \boxed{ \nabla^2\mathbf{B}_\text{NMR} = \frac{\mu_0}{3}\nabla^2\mathbf{M} -\mu_0\nabla\left(\nabla\cdot \mathbf{M}\right),} \end{align} where \(\nabla^2=\nabla\cdot\nabla\) denotes the Laplace operator.

Since \(\chi \le1\), we have the following approximation

\[\begin{align} \mathbf{M} = & \frac{1}{\mu_0} \frac{\chi}{1+\chi} \left(\mathbf{B_0} + \delta \mathbf{B}\right) \\ \simeq & \frac{\chi}{\mu_0}\mathbf{B_0} \end{align}\]

Here we have used the fact that \(\delta B\ll B_0\).

scalar susceptibility

Assuming scalar susceptibility, the signal equation can be simplified into \begin{align} \boxed{ \nabla^2 f_T= \left(\frac{\nabla^2}{3}-\partial_z^2\right)\chi \equiv \Box \chi,}
\end{align} where \(f_T= (B_\text{NMR}-B_0)/B_0\) is the relative difference field (RDF).

If \(\chi\) is known, the relative field deviation \(f_T\) is given by \begin{align} f_T = -\frac{1}{4\pi} \int \frac{\Box \chi}{|\mathbf{r}-\mathbf{r’}|} d\mathbf{r’} \end{align} since the Green’s function of the Laplacian operator over the full space is \(1/r\), i.e., \begin{align} \nabla^2 \frac{1}{r} = -4\pi \delta(\mathbf r). \end{align}

Using integration by parts, we get \begin{align} f_T(\mathbf r) = \lim_{\epsilon\rightarrow 0} \int_{|\mathbf{r}-\mathbf{r’}|>\epsilon} d(\mathbf{r}-\mathbf{r’}) \chi(\mathbf{r’})d\mathbf {r’}\equiv d*\chi. \end{align} where \(*\) denotes convolution and the dipole field along z direction is \begin{align} d(\mathbf r) = \frac{1}{4\pi}\frac{3\cos^2\theta -1}{r^3}. \end{align} Note \(d(0)\) needs to be excluded.

So far the equations are all written in \(\mathbb R^3\), in reality, \(f_T\) is only known in a localized region \(\Omega\). Thus the measured \(f_T\) can be separated into two parts, \(f_B\) caused by sources \(\chi_B\) outside \(\Omega\) and \(f_L\) caused by sources \(\chi_L\) inside \(\Omega\), i.e., \begin{align} f_T(\mathbf r) = f_B(\mathbf r) + f_L(\mathbf r), \mathbf r \in \Omega \end{align} and \(f_L = d* \chi_L\) and \(f_B = d*\chi_B\) where \(\chi_B\) is defined on \(\bar\Omega\) and \(\chi_L\) is defined on \(\Omega\).

Since \(\nabla^2 d = 0\) on \(\mathbb{R}^3\backslash\{0\}\) and \(f_B\) , we have \begin{align} \nabla^2 f_B = 0. \end{align} If the boundary condition of \(f_B\) is known on \(\partial\Omega\), one can either solve the boundary value problem of the Laplace’s equation or match boundary values with \(f_B=d*\chi_B\) to remove background field.

inverse problem: regularization

The QSM problem is to solve \(\chi\) given \(f_T\) according to the signal equation. The difficulty of this inverse problem is two folds. Firstly, the boundary values of \(\chi\) are not available. Secondly, the differential operator \(\Box\) has a non-trivial null space.

In practice, regularized minimization is used to solve this inverse problem \begin{align} \min_{\chi} |f_L-d*\chi| + \lambda |reg|. \end{align}

In morphology enabled dipole inversion (MEDI), the minimization problem is \begin{align} \min_{\chi} |f_L-d*\chi|_2 + \lambda |\nabla \chi/\nabla I|_1. \end{align} Here the regularization term enforces the insight that a boundary is unlikely to exist in the susceptibility map if there is no boundary at the at location in the reference intensity image. \(1/\nabla I\) is implemented as a binary mask.

inverse problem: solving PDE

The QSM signal equation takes the form of 2D wave equation with source if the z direction is treated as time. It can be solved directly if intial-boundary conditions are available. In practice, one can make judicious guess based on prior informations or simply use homogeneous initial-boundary conditions. On the other hand, \(f_T\) is derived from experimental measurements thus the source term contains error due to noise. It is then important to study the error induced by the errors in the source term as well as errors in the initial-boundary conditions.

The error \(\epsilon(x,t)\) satisfies the following equation

\[\begin{align} \begin{cases} \left(\partial_t^2+D\right)\epsilon(x,t) = \epsilon_f(x,t) \\ \epsilon(x,0) = \epsilon_1(x) \\ \epsilon'(x,0) = \epsilon_2(x) \end{cases} \end{align}\]

where \(D\) is a spatial differential operator, \(\epsilon_{1,2}(x)\) denote the error in the initial-boundary conditions. For simplicity, we assume the error in the source term \(\epsilon_f\) takes the following form \begin{align} \epsilon_f(x,t) = \sum_k A_k \delta(x-x_k)\delta(t-t_k). \end{align}

To solve it, we utilize the spectral decomposition of \(D\) \begin{align} D v_i(x) = \lambda_i v_i(x). \end{align} Here \(v_i(x)\) form a set of spatially orthonormal basis functions. It is easy to show that \begin{align} \epsilon(x,t) = \sum_iv_i(x)\left[\epsilon_{1i}\cos\sqrt{\lambda_i}t +\frac{\epsilon_{2i}}{\sqrt{\lambda_i}}\sin\sqrt{\lambda_i}t + \epsilon_{3i}(t) \right], \end{align} where \(\epsilon_{1i}\equiv\left<v_i,\epsilon_1\right>\) is the projection of \(\epsilon_1(x)\) on the \(i\)‘th basis. \(\epsilon_{2i}\) is defined accordingly. \(\epsilon_{3i}\) has the form \begin{align} \epsilon_{3i} = \sum_k\frac{v_i(x_k)}{\sqrt{\lambda_i}} \sin\sqrt{\lambda_i}(t-t_k)\theta(t-t_k), \end{align} where \(\theta(t)\) is the unit step function.

generalized Lorentz correction for white matter regions

In white matter regions, the Lorentz sphere correction does not hold and we use a different model \begin{align} f = d*\chi + \ell \chi, \end{align} where \begin{align} \ell = \left(\frac{1}{3} - \frac{\sin^2\theta}{2}\right)M_\text{WM} \end{align} is a diagonal matrix and \(\theta\) is the angle between \(B_0\) direction and the fiber direction. The white matter mask \(M_\text{WM}\) can be generated by thresholding the FA map.

The energy function takes the form \(\|w(Q\chi-f)\|^2\). Here \(w\) is a weighting function, \(Q\equiv F^{-1}DF+\ell\) with \(F\) is the Fourier transform matrix, and \(D\equiv Fd\) is the dipole function in Fourier space. It gives rise to the linear equation \begin{align} Qw^2Q\chi = Qw^2f. \end{align}