The-Vinh Tran-LuuUniversity of Maryland, College Park, MD, 20742, USAMark-Alexander HennUniversity of Maryland, College Park, MD, 20742, USANational Institute of Standards and Technology (NIST), Gaithersburg, MD, 20895, USAKlaus N. QuelhasNational Institute of Standards and Technology (NIST), Gaithersburg, MD, 20895, USASolomon I. WoodsNational Institute of Standards and Technology (NIST), Gaithersburg, MD, 20895, USA

(July 30, 2024)

###### Abstract

Magnetic particle imaging (MPI) is an emerging imaging technique with many applications and a very active field of research. This app provides users with the opportunity to develop some intuition about the inner workings of MPI as it is being researched through NIST’s Thermal MagIC project in an interactive and fun way. Users can vary different experimental and post-processing parameters to see how the image quality and particle reconstruction changes for different measurement conditions.

## 1 Introduction

Magnetic particle imaging (MPI) was introduced almost twenty years ago as a way for remotely detecting magnetic nanoparticle (MNP) tracers, with several applications in biomedical imaging and diagnosis, as well as materials research [8, 4, 2]. MPI relies on the non-linear magnetization response of MNPs when exposed to time-varying magnetic fields. By generating a moving field-free point (FFP) which saturates all particles but the ones near the no-field region we can scan a specified field of view and measure the magnetization response of the particles in the region of the FFP to a time-varying excitation field.

However, the time signal of the varying magnetic response alone is not enough to deduce the spatial distribution of the MNPs. Obtaining this information requires two additional steps, in which we first relate the time signal to its corresponding spatial position to get an image, and then deblur this image based on the knowledge of the experimental setup and its point spread function (PSF). This process in mathematical terms is also known as solving an inverse problem and amounts to carefully balancing between a good fit to the observed data and a realistic reconstruction.

This web application intends to give the user a tool to see for themselves how different experimental setups can influence the quality and reliability of Magnetic Particle Imaging and can help to illustrate the challenges researchers at NIST face when trying to improve this novel measurement technology[1, 3, 6]. It employs a simplified version of NIST’s MPI library for Python^{1}^{1}1Reference is made to commercial products to adequately specify the experimental procedures involved. Such identification does not imply recommendation or endorsem*nt by the National Instituteof Standards and Technology, nor does it imply that these products are the best for the purposespecified., and runs using Python Flask, HTML, and CSS, as well as matplotlib and other dependencies. The project code is avaliable for viewing and download on GitHub. You can run the project by cloning the repository, installing the dependencies in the requirements.txt text file, and running the command ”python app.py” in the project directory. Alternatively, the user can access a web-based version of the app via this link.

Similar to MPI, the app consists of two major parts: the signal generation, and the reconstruction of the particle distribution from the generated signal.

### 1.1 Signal Generation

The signal measured in an MPI setup can be modeled using:

$s(t)=-B_{s}\frac{d}{dt}\int_{\Omega}\rho(x)m\mathcal{M}_{\textsf{eff}}[H_{%\textsf{app}}(x,t)]dx,$ | (1) |

with $\rho$ the particle distribution, $m$ is the magnetic moment of a single particle given by the product of the saturation magnetization $M_{s}$ in [A/m] and the volume of the core of a single particle, that itself depends on the particle’s diameter $d$, $B_{s}$ is the coil sensitivity given in [T/A], and $\mathcal{M}_{\textsf{eff}}$ is a function describing the magnetization response to a time-varying applied field $H_{\textsf{app}}$.

A commonly applied model for this magnetization behaviour is:

$\displaystyle\mathcal{M}_{\textsf{eff}}[H_{\textsf{app}}(x,t)]$ | $\displaystyle=$ | $\displaystyle\mathcal{L}[\kappa\|H_{\textsf{app}}(x,t)\|]\frac{H_{\textsf{app}%}(x,t)}{\|H_{\textsf{app}}(x,t)\|},$ | (2) |

with the so-called Langevin function

$\displaystyle\mathcal{L}[\kappa\|H_{\textsf{app}}(x,t)\|]$ | $\displaystyle=$ | $\displaystyle\coth[\kappa\|H_{\textsf{app}}(x,t)\|]-\frac{1}{\kappa\|H_{%\textsf{app}}(x,t)\|},$ | (3) |

and $\kappa=\frac{\mu_{0}m}{k_{B}T}$, $k_{B}$ is the Boltzmann constant and $T$ is the temperature in $K$. MPI utilizes the non-linear magnetization response by creating a field-free point (FFP) that is sensitive to high-frequency changes of the applied field. This is done by combining a gradient field with a gradient $G$ matrix, and a so-called drive field $H_{d}$, such that the position of the FFP is given by

$r(t)=G^{-1}H_{d}(t).$ | (4) |

In the present case, the drive field is defined via

$H_{d}=\begin{bmatrix}H_{0}\sin(2\pi f_{0}t)\\H_{1}\sin(2\pi f_{1}t)\end{bmatrix},$ | (5) |

which leads to a so-called Lissajous trajectory for the FFP.

From Eq. (4) it follows that the combination of the gradient field and the drive-field amplitudes determines the maximum range of the trajectory $r$, and hence the dimensions of the field of view (FOV). In the present case, we fix the dimensions of the FOV, such that the drive field amplitudes $H_{0}$ and $H_{1}$ become dependent variables. Using the information about the FFP trajectory, we can correlate the time signal from Eq. (1) with its spatial position and obtain what sometimes is called the raw image: a blurred version of the original particle distribution, with the degree of blurring determined by the experimental parameters, and ultimately the point spread function of the system.

### 1.2 Image Reconstruction

The MPI system is a linear-shift invariant system [5]. Hence, the $n$-dimensional signal vector $y$ stemming from an arbitrary distribution of particles $\rho$ can be understood as a weighted sum of signals from $\delta$-samples placed at different positions within the FOV that has been discretized into $m$ voxels. The weights are proportional to the number of particles at that position, and can be determined from a simple matrix equation:

$S\rho=y+\epsilon,~{}\mathrm{with}~{}S\in\mathbb{R}^{m\times n},\rho\in\mathbb{%R}^{m},y\in\mathbb{R}^{n}.$ | (6) |

The matrix $S$ is called the system matrix, its columns correspond to the measured $n$-dimensional signal from a $\delta$-sample put at the $m$ different positions within the FOV, the term $\epsilon$ denotes possible noise in the measurement. In general $m<n$, which means we are faced with an overdetermined system, so the inverse of the matrix $S$ can not be defined in a straight-forward manner, and instead of solving for $\rho$ directly we need to consider the optimization problem [7]:

$\hat{\rho}=\underset{\rho\in\mathbb{R}^{m}}{\mathrm{argmin}}\|S\rho-(y+%\epsilon)\|.$ | (7) |

Since this system of linear equations is overdetermined, there can be more than one solution that fits the measurement data well. It is however reasonable to expect the solution to this equation to be somewhat regular or smooth, or in other words to have a small $l^{2}$-norm, leading us to add a so-called regularization term, given by $\lambda_{r}\|\rho\|$ to our cost function. This leads to the regularized version of our problem:

$\hat{\rho}=\underset{\rho\in\mathbb{R}^{m}}{\mathrm{argmin}}\left\{\|S\rho-(y+%\epsilon)\|+\lambda_{r}\|\rho\|\right\}.$ | (8) |

By varying the value of $\lambda_{r}$ we can balance between the goodness-of-fit to the signal and the desired smoothness of the solution.

## 2 User Guide

Our web-based app employs the models and techniques introduced in the previous section to provide an interactive environment that can help to develop a better understanding of the MPI framework. This process again consists of two separate steps: the signal generation and the reconstruction.

### 2.1 Signal Generation in App

The first interactive screen is the simulation screen, in which the raw time signal data can be generated. Once the user hits the Run Code! button the app starts calculating the necessary data and puts out four plots. These consist of the ground truth particle distribution and the magnetization curve, see Fig. (2), and the FFP trajectory and the raw signal, see Fig. (3). Noise is added to the generated data, via the term $\epsilon$ in Eq. (6). It is drawn from a normal distribution with zero mean and a standard deviation that is proportional to the maximum of the generated signal, hence $\epsilon\sim\mathcal{N}(0,\sigma^{2})$, with $\sigma=\kappa\cdot\max(y)$, such that a noise level of $0.01$ corresponds to a standard deviation of 1% of the maximum of the measured signal.

Parameter | Symbol | Default value |
---|---|---|

Gradient field in $x$-direction | $g_{00}$ | 8 T/m |

Drive coil frequency in $x$-direction | $f_{0}$ | 25.5 kHz |

Drive coil frequency in $y$-direction | $f_{1}$ | 25.25 kHz |

Particle diameter | $d$ | 20 nm |

Saturation magnetization | $M_{s}$ | $4.5\cdot 10^{5}$ A/m |

Coil sensitivity | $B_{s}$ | $8\cdot 10^{-4}$ T/A |

Regularization parameter exponent | $\lambda$ | -6 |

Noise level | $\kappa$ | $10^{-3}$ |

Note, that the ground truth particle distribution is not affected by varying the experimental parameters, as are the $x$ and $y$ dimensions of the FOV. The M vs. H curve changes with $M_{s}$, $d$, and $g_{00}$ (see Table 1) such that an increase in $M_{s}$ or $d$ leads to a steeper shape, and a change in only $g_{00}$ only changes the range of $H$ over which we plot the M vs. H curve. It is important to note that the user only provides a single value $g_{00}$ for the gradient matrix G, from we which construct the gradient matrix as:

$G=\begin{bmatrix}g_{00}&0\\0&-\frac{1}{2}g_{00}\end{bmatrix}$ |

In our setting, the FFP trajectory is only a function of $f_{0}$ and $f_{1}$, since we fix the FOV and change the drive field amplitudes accordingly, while the raw signal depends on all mentioned parameters, and the noise level. Note, that changing the value of the coil sensitivity $B_{s}$ only scales the generated signal.

By clicking on the Next (Reconstruction) button, the user advances to the next screen which provides the reconstructed particle distribution along with some additional information.

### 2.2 Image Reconstruction in App

Using the data generated in the previous screen, we determine the particle distribution from the raw signal data, and are presented with five plots. The first two in Fig. (5) show the time signal related to its spatial position, also called the raw image, and the point spread function (PSF) which gives us an idea about the blurring that stems from the selected experimental parameters. In general a larger value for $g_{00}$ leads to less dispersed PSF, and hence a better resolution in the raw image. The second set of plots, Fig. (6), show the reconstructed particle distribution in terms of profiles along the two FOV axes, and a 2D representation of the reconstruction.

When you hit the Run Again! button the reconstruction is rerun. This is helpful to investigate how the regularization parameter $\lambda$ or various experimental parameters change the reconstruction. Note, that the effective regularization parameter $\lambda_{r}$ in Eq. (8) is calculated as $\lambda_{r}=10^{\lambda}$.

## References

- [1]T.Q. Bui, M.-A. Henn, W.L. Tew, M.A. Catterton, and S.I. Woods.Harmonic dependence of thermal magnetic particle imaging.Scientific Reports, 13(1):15762, 2023.
- [2]B.Gleich and J.Weizenecker.Tomographic imaging using the nonlinear response of magnetic particles.Nature, 435(7046):1214–1217, 2005.
- [3]M.-A. Henn, K.N. Quelhas, T.Q. Bui, and S.I. Woods.Improving model-based MPI image reconstructions: Baseline recovery, receive coil sensitivity, relaxation and uncertainty estimation.International Journal on Magnetic Particle Imaging, 8(1), 2022.
- [4]T.Knopp and T.M. Buzug.Magnetic Particle Imaging: An Introduction to Imaging Principles and Scanner Instrumentation.Springer, Berlin/Heidelberg, 2012.
- [5]K.Lu, P.W. Goodwill, E.U. Saritas, B.Zheng, and S.M. Conolly.Linearity and shift invariance for quantitative magnetic particle imaging.IEEE transactions on medical imaging, 32(9):1565–1575, 2013.
- [6]K.N. Quelhas, M.-A. Henn, R.C. Farias, W.L. Tew, and S.I. Woods.GPU-accelerated parallel image reconstruction strategies for magnetic particle imaging.Physics in Medicine & Biology, 69(13), 2024.
- [7]A.Tarantola.Inverse problem theory and methods for model parameter estimation.SIAM, 2005.
- [8]J.Weizenecker, J.Borgert, and B.Gleich.A simulation study on the resolution and sensitivity of magnetic particle imaging.Physics in Medicine & Biology, 52(21):6363, 2007.