Towards Accurate Modeling of the Multidimensional Magnetic Particle Imaging Physics

The image reconstruction problem of the tomographic imaging technique magnetic particle imaging (MPI) requires the solution of a linear inverse problem. One prerequisite for this task is that the imaging operator that describes the mapping between the tomographic image and the measured signal is accurately known. For 2D and 3D excitation patterns, it is common to measure the system matrix in a calibration procedure, that is both, very time consuming and adds noise to the operator. The need for measuring the system matrix is due to the lack of an accurate physical model that is capable of describing the nanoparticles' magnetization behavior. Within this work we introduce a physical model that is based on N\'{e}el rotation for large particle ensembles and we find model parameters that describe measured 2D MPI data with much higher precision than state of the art MPI models. With phantom experiments we show that the simulated system matrix can be used for image reconstruction and reduces artifacts due to model-mismatch considerably.


A. Introduction
Tomographic imaging techniques play a major role in medical diagnostics and treatment. In 2005, a new method called magnetic particle imaging (MPI) was introduced [3], which allows to image magnetic nanoparticles (MNPs) with high sensitivity and high spatio-temporal resolution. In contrast to MRI, the contrast in MPI is positive and the reconstructed images are quantitative. With these properties MPI is an interesting candidate for various medical applications ranging from vascular imaging [23] to targeted imaging [32].
Signal encoding in MPI is achieved by exciting the MNP with one or several homogeneous excitation fields with a frequency in the kHz region [11]. The MNP respond with a change of their magnetization, which is measured with one or multiple receive coils. Spatial encoding is achieved by superimposing a static gradient field leading to a spatially dependent magnetization response of the MNP to the excitation field. In ferrofluids the change of the particle magnetization is enabled by two different dynamic processes [6]: The Brownian rotation describes the mechanical alignment of the entire particle with a change of the magnetic field whereas Néel rotation describes the alignment of the particle's inner magnetic moment. Note that due to the continuous excitation the system does not relax in MPI. These dynamics in the context of MPI were initially studied by Weizenecker et al. [26] in terms of stochastic ordinary differential equations. A first attempt to Brownian rotation via the Fokker-Planck equation was presented by Yoshida and Enpuku for 1D excitation patterns [28] followed by further works in this direction [27,29,31,20,19]. Different MNP are often categorized by the dominating mechanism. For example, Resovist's magnetization behavior is mainly determined by Néel rotation, whereas a frequency-dependent (with respect to the applied field) influence of Brownian rotation can be observed [14]. Furthermore, the orientation of the particles' easy axis significantly influences the magnetization behavior due to the Néel rotation [30]. The still incomplete nature of models for imaging in MPI has motivated an increasing number of studies in the context of Brownian and Néel rotations for MPI excitation patterns [27,15,1,2,21,5].
In order to reconstruct an image of the MNP concentration, the physical processes during the MPI experiment need to be inverted. Mathematically, this involves the solution of an ill-posed linear system of equations [25]. One prerequisite for image reconstruction is that the system matrix is accurately known. Since its very first introduction in [3], the system matrix is measured in a data-driven calibration procedure  Figure 1: Imperfections of the equilibrium model compared to measured data. Absolute of two selected frequency components (denoted by mixing factors) are shown in the first and second row. The measured system matrix using the unmodified particle solution is shown in the first column while the simulated matrix based on the equilibrium model (model A) is shown in the second column. The sampling pattern is shown in the last column.
by determining the response of the system to a small delta sample that is shifted through the imaging volume. The main advantage of this method is that it takes all physical effects of the imaging process into account but it has also two major drawbacks: The measured system matrix includes noise and the calibration procedure is very time consuming lasting up to several days [22]. These drawbacks motivated the development of methods determining the MPI system matrix based on physical models. While the initial model-based results in [13] and [10] were promising, model-based reconstruction based on the used equilibrium model denoting the simplified model which assumes infinitely fast relaxation leads to worse image quality than the data driven approach [8]. This is why the latter is still the method of choice in almost all publications since 2010 that use multidimensional excitation patterns. For 1D excitation, which includes 2D and 3D Cartesian like sampling patterns [9],model-based reconstruction has been established in various works since it enables simple time-signal based reconstruction [4]. Within this work we investigate the question why the simple equilibrium model fails to accurately describe the MPI system matrix for 2D Lissajous type excitation patterns. As a motivating example we consider a system matrix measured along a 2D Lissajous trajectory (frequency ratio 16/17, details are outlined in Section F). Fig. 1 shows two selected rows of the 2D MPI system matrix and compares the measurement with the equilibrium-based simulation. One can see that the simulation matches quite well the basic structure of the matrix rows, which are consisting of wave pattern that share high similarity to 2D tensor products of Chebyshev polynomials [17]. However, when taking a closer look one can identify various differences that lead to larger numerical deviations. In particular we want to highlight two qualitative effects. First, the wave hills are merging in outer regions of the field-of-view (FOV), i.e. in regions where the simulation has a zero-crossing there is a non-zero value in the measured one. Second, one can see a tilting of the outer wave structures in outer regions of the FoV. The goal of this work is to simulate these effects using a physical model that takes Néel rotation of the particles' magnetic moments into account.

B. Methods
We consider a general MPI imaging experiment where the particles are located in the MPI scanner and the change of the particle magnetization is measured using induction coils. Then the signal induced in a receive coil with sensitivity p : where c : Ω → R + 0 is the concentration of the magnetic nanoparticles,m : R 3 × [0, T ] → R 3 is the mean magnetic moment of particles and Ω ⊂ R 3 is the imaging volume. The mean magnetic momentm(x, t) depends on the applied magnetic field H app : R 3 × [0, T ] → R 3 , which is usually a T -periodic function with period T along the time dimension. The induced signal is then convolved with the T -periodic analog filter a : R → R to remove the main contributions of the direct feedthrough yielding the signal v = a * ṽ. Therefore, v is T -periodic as well and can be expanded into a Fourier series with coefficientŝ This formulation is commonly used, since the signal at the excitation frequencies is blocked using an analog band-stop filter prior to the signal digitization. With where s k : R 3 → C is the system function that we already discussed in Fig. 1. When sampling the system function s k at discrete positions x l , l = 1, . . . , N , and with a certain sampling rate in time one obtains the MPI system matrix S = (s k (x l )) k=0,...,K;l=1,...,N . In this work we will consider both modeled system matrices S Mod as well as measures system matrices S Cal that are obtained by moving a small delta sample through the FoV. In case of multiple receive channels, the corresponding system matrices can be stacked to arrive at a joint linear system of equation; see for example [7] for a formal definition.

C. Particle Models
Since H app and p can be easily modeled using the Biot-Savart law and a can either measured or fitted [10], we focus on modeling the particles' mean magnetic moment behavior. We first consider the equilibrium model that can be derived under the assumption that the particles are always in thermal equilibrium in which casē m directly follows the applied field H app . We refer to it as model A which is given bȳ where L β : R → R is given in terms of the Langevin function by the following: kBTB mainly depend on the particle core diameter D. Since we have already seen that the equilibrium model is not accurate enough for modeling the physics in an MPI experiments we next consider a more appropriate model that takes dynamic relaxation effects into account. There are basically two different approaches, one can either consider the problem on a micromagnetic level for individual particles and solve the Langevin equation corresponding to the Landau-Lifschitz-Gilbert equation for a sufficiently large number of particles to obtain a reasonable estimate for the mean. Alternatively, one can take a comprehensive view and solve the Fokker-Planck equation for a probability distribution representing an entire ensemble of nanoparticles in terms of a parabolic partial differential equation. The latter approach benefits from efficient and more accurate solutions when aiming for a mean computation [24] as determining the mean from the Langevin equation requires a large number of particle simulations.
The second case which we refer to as model B is based on the latter approach which determines the mean magnetic moment via the probability density function f : which is the solution to the corresponding Fokker-Planck equation where S 2 is the surface of the sphere in R 3 . The mean is then given bym where f is the solution to the following specific case of a convection-diffusion equation on the sphere where τ > 0 is the relaxation time constant and the (velocity) field a : where p i ≥ 0, i = 1, . . . , 4, and n ∈ S 2 is the easy axis of the particles (or alternatively n : Ω → S 2 and p 3 , p 4 : Ω → R + 0 ). Differentiation in terms of gradient ∇ S 2 and divergence div S 2 is considered with respect to the surface S 2 . (·, ·) denotes the Euclidean scalar product of R 3 . The specific structure of the Fokker-Planck equation is due to the independent consideration of Brownian and Néel rotation model as it was observed for common tracer that one meachism can be dominant [14]. For a detailed discussion of the relationship between Fokker-Planck equation and Langevin equations for Brownian and Néel rotation we refer to the survey [6]. A pure Néel rotation including anisotropy is given by p 1 =γµ 0 , p 2 =γαµ 0 , p 3 = 2γ Kanis MS , p 4 = 2γ Kanis MS , and τ = VCMS 2kBTBγα (γ = γ 1+α 2 ). We note that the parabolic PDE in (4) has no dependence on derivatives with respect to the spatial variable x. It can thus be considered as parametric with respect to x.
In the following we use the Néel rotation model as it includes the particle anisotropy. Here we distinguish the following three cases to identify the relevant modeling aspects: B1 Néel rotation model without anisotropy, i.e., p 3 = p 4 = 0 in (5). B2 Néel rotation model including anisotropy, i.e., it includes all summands in (5) for a given easy axis n ∈ S 2 . B3 Néel rotation model including a space-dependent anisotropy which can be motivated by the local structure of the magnetic field, i.e., n : Ω → S 2 . In particular in the corners of the FoV the magnetic field vector is mainly oriented within one quadrant. During multiple repetitions of the excitation sequence this can lead to a preferred orientation of the easy axis of the nanoparticles in a ferrofluid, which can highly influence the magnetization behavior as reported in [30] Following the approach in [24] the numerical solution of (4) is obtained as follows: We consider an initial value f (m, x, 0) = f 0 (m, x) with f 0 (·, x) ≥ 0 and f 0 (·, x) L 1 (S 2 ) = 1 for any x ∈ Ω. A semi-discrete Cauchy problem is obtained via discretization in non-normalized spherical harmonics {Y m l } l,m (with respect to m while considering a rotated problem using a rotation matrix R ∈ R 3×3 such that Rn = e 3 ) which yields a first order ordinary differential equation system for any x ∈ Ω. The initial value problem is then solved on a time interval covering multiple periods of the excitation pattern while using a uniform initial value (for any x ∈ Ω). The numerical solution is obtained via a variable order, variable step solver using Matlab's builtin function ode15s. The last period of the solution is then used for further investigations. The mean magnetic moment can then be determined directly via linear combinations of the coefficients corresponding to the spherical harmonics Y 0 1 , Y 1 1 , and Y −1 1 , which is subsequently rotated back.

D. Distance Measure
For quantifying the accuracy of a model we express the frequency index k in terms of the mixing orders m x , m y [18] fulfilling k(m x , m y ) = m x N Dens + m y (N Dens + 1) Here, N Dens is the parameter that controls the density of the Lissajous sampling pattern and for each k the factors m x and m y with the smallest absolute value fulfilling (6) are considered. Now with S Cal k(mx,my) being the measured and S Mod k(mx,my) being the modeled system matrix row we consider the the normalized root mean square error (NRMSE) In addition to this per-row metric we consider the mean NRMSE over the m x = 0, . . . , M x and m y = 0, . . . , M y where M x , M y ∈ N are upper bounds for the mixing order.

E. Image Reconstruction
For image reconstruction we use the standard framework [12] that applies the iterative Kaczmarz algorithm augmented by Tikhonov regularization and a positivity constraint. We apply standard SNR thresholding where matrix rows with low SNR (¡ 2.0) are removed prior to reconstruction. In addition we apply a second way of matrix row filtering. The idea is to select only those matrix rows of the modeled system matrices, which deviate only marginal from the measured one. Let Θ ∈ [0, 1] be a predefined threshold. Then, we choose matrix rows k where ε(S Cal k , S Mod k ) < Θ. The parameter Θ can now be used to trade of between inaccuracies due to the model/measurement mismatch and a loss off spatial resolution that happens when only few system matrix rows are selected for reconstruction.

F. Experimental Setup
Experimental data was measured with a preclinical MPI scanner (Bruker, Ettlingen) using field-free-point (FFP) excitation patterns. In all experiments a 2D cosine excitation (N Dens = 16) within the xy-plane of the scanner was applied with excitation-field amplitudes A x = A y = 12 mTµ −1 0 and 1 Tm −1 µ −1 0 selection field gradient within the xy-plane. The resulting area covered by the FFP is of size 24 × 24 mm 2 .
The system matrix was measured by shifting a 1 mm 3 delta sample of Perimag (Micromod, Rostock) to 30 × 30 positions covering a FoV of size 30 × 30 mm 2 (1 mm step size). A particle phantom consisting of three rods 1.0 mm in diameter is prepared using the same tracer of the same concentration as the delta sample. The phantom is shown in Fig. 6.

G. Model Parameter Selection
The considered models each have several parameters that are a-priori unknown. Since the simulation of the models requires a significant time, it is challenging to apply for example gradient-based optimization techniques to the continuous parameter space. Therefore, we calculate system matrices for a finite number of settings and perform a subsequent discrete optimization.
Model A as well as models B1-B3 require as input the particle core diameter D. We simulated the equilibrium model for D ∈ {20 nm, 25 nm, 30 nm} and found in initial tests that the system matrices with D = 20 nm fitted best in all model cases and are therefore considered for further investigation. For the Néel model B2 we consider a spatially homogeneous easy axis distribution n(θ) = (cos(θ), sin(θ)) T and orientations θ ∈ {0 • , 45 • , 90 • , 135 • }. In this case the anisotropy constant is chosen as K anis = 625 J/m 3 .
For the inhomogeneous case B3 we consider an easy axis n(x) = x |x| and K anis (x) = g Kanis |H S (x)| where we chose g Kanis ∈ {625/h, 1250/h, 2500/h, 5000/h}J/m 3 with h = max x∈Ω |H S (x)|. Thus, the easy axis is parallel to the selection field vector and the effective anisotropy increases when increasing the distance to the FFP of the selection field.

H. Results
First, we consider the Néel model B2 with spatially homogeneous anisotropy distribution and homogeneous easy axis. Selected rows of simulated system matrices are shown in Fig. 2. One can see that the angles 0 • and 90 • do not change the shape of the patterns compared to the equilibrium model. However, for 45 • and 135 • one can see the desired effect that the patterns are tilted in one direction, which is the orientation of the easy axis. In addition one can also see the merge of the wave hills that was observed in the measured system matrices. However, in direct comparison with the measured system matrix, one can see that the direction of the tilt is homogeneous while the measurement shows a spatially dependent tilt. One attempt could be to assume that there is a homogeneous distribution of the easy axis in each voxel, which can be simulated by taking the mean over all simulated angles. But as is shown in Fig. 2, applying the mean leads to a similar pattern as observed for the equilibrium model A (see Fig. 1).
The results from the spatially homogeneous anisotropy distribution motivate the usage of an spacedependent anisotropy distribution. In order to make the tilting orientation-dependent, we let the easy axis   4). Column 5 shows the mean over all easy axes while column 6 shows the measured system matrix for reference.  Figure 3: Selected frequency components of modeled system matrices using the Néel model B3 with anisotropy gradients ranging from 5000/h J/m 3 to 625/h J/m 3 (columns 1-4;left to right). Column 5 shows the modeled system matrix using the Néel model B1 without anisotropy while column 6 shows the measured system matrix for reference. Below, the mean NRMSE between the model and the measurement is shown.  Model A Model B3 Figure 6: Reconstruction results of measured phantom data using three different system matrices: The measured system matrix, the matrix based on the equilibrium model A, and the system matrix obtained from the Néel model B3 for g 3 Kanis = 1250/h J/m 3 . Each reconstruction uses the same regularization parameter λ = 10 −7 and the threshold Θ is varied between 0.12 and 0.3. introduce the two effects observed in the measurements. First a tilting effect and second a merging of the wave hills. With increasing anisotropy gradient, the effect gets stronger and starts already at the center. For small gradients, the effect is only visible in outer regions. For g 3 Kanis = 1250/h J/m 3 we find the highest similarity compared to the measurement, which is also supported by the mean NRMSE shown in Fig. 3.
After choosing and optimizing the model, we next make a more detailed comparison of model and measurement. Fig. 4 shows various system matrix rows of the optimized Néel model B3 in comparison with the equilibrium model and the calibration for various mixing factors and the x receive channel. One can see that the optimized Néel model much better describes the measurement than the equilibrium model. To support this observation, the figure also shows difference images. Just for m x = 0 and m x = 1, 2 the pattern still looks slightly different and the error is similar for both models. To quantify the differences we calculate the NRMSE for m x = 0, . . . , 9 and m y = 0, . . . , 9 which are shown in Fig. 5. The data confirms that the Néel model has lower error than the equilibrium model for most of the mixing factors. In particular for m x > m y the error is much lower. For m x m y , however, the error increases and there is no clear advantage of the Néel model.
Reconstruction results for the phantom data are shown in Fig. 6. One can see that the Néel model achieves a similar image quality as the calibrated system matrix. As expected the results are closer to each other if the threshold Θ is decreased. For large thresholds, the Néel model B3 results in artifacts at the boundaries of the FoV. The equilibrium model A shows a much worse image quality, which underlines, why it is usually not used in practice.

I. Discussion
The results show that the physics of a multi-dimensional MPI experiment can be modeled by solving the Fokker-Planck equation considering Néel rotation with a spatially inhomogeneous anisotropy. Since the simulation of the Néel model is computational expensive we did not fully optimize the set of unknown parameters yet. One can expect that the accuracy of the model increases even further when considering a particle size distribution and optimizing the probability density function of the particles. Instead of the ad-hoc choice of a linear increasing gradient anisotropy it seems also promising to optimize the shape of the anisotropy field. The next step for a more precise model could be to consider a coupled Brownian and Néel rotation model as for example has been derived in [24].
The strengths of using a model compared to using a calibration are manifold. We see the prime advantage that it allows to decouple the magnetic field sequence from the particle behavior. This enables to use the same model and evaluate it under different field conditions (different excitation-field amplitudes and gradient strength). The need for such an approach is even more pressing for multi-patch sequences where a multitude of similar but slightly perturbed system matrices need to be acquired since the fields show imperfections in space [22].
Another application is multi-contrast MPI [16]. Here, one requires system matrices for different states of the particles (e.g. different viscosity or temperature). Using an accurate model, it is possible to determine the required system matrices just by modifying the parameters in the model. Furthermore, the problem of multi-contrast MPI can be formulated in terms of a joint parameter identification problem.