Modified Jiles–Atherton Model for Dynamic Magnetization in X-Space Magnetic Particle Imaging

Objective: Magnetic particle imaging (MPI) is a promising medical modality that can image superparamagnetic iron-oxide nanoparticle (SPIO) concentration distributions safely and with high sensitivity. In the x-space reconstruction algorithm, the Langevin function is inaccurate in modeling the dynamic magnetization of SPIOs. This problem prevents the x-space algorithm from achieving a high spatial resolution reconstruction. Methods: We propose a more accurate model to describe the dynamic magnetization of SPIOs, named the modified Jiles–Atherton (MJA) model, and apply it to the x-space algorithm to improve the image resolution. Considering the relaxation effect of SPIOs, the MJA model generates the magnetization curve via an ordinary differential equation. Three modifications are also introduced to further improve its accuracy and robustness. Results: In magnetic particle spectrometry experiments, the MJA model shows higher accuracy than the Langevin and Debye models under various test conditions. The average root-mean-square error is 0.055, 83% and 58% lower than the Langevin and Debye models, respectively. In MPI reconstruction experiments, the MJA x-space improves the spatial resolution by 64% and 48% compared to the x-space and Debye x-space methods, respectively. Conclusion: The MJA model shows high accuracy and robustness in modeling the dynamic magnetization behavior of SPIOs. By integrating the MJA model into the x-space algorithm, the spatial resolution of MPI technology was improved. Significance: Improving the spatial resolution by using the MJA model helps MPI have a better performance in medical fields, including cardiovascular imaging.


I. INTRODUCTION
M AGNETIC particle imaging (MPI) is a promising medical modality that can image the concentration distribution of superparamagnetic iron-oxide nanoparticles (SPIOs) safely and with high sensitivity [1]. Thus, MPI has considerable application potential in many fields, including cancer detection [2], cell therapies [3], cardiovascular imaging [4], [5], and stem-cell tracking [6]. This technology exploits a gradient field to create a field-free point (FFP) or field-free line (FFL) [1], [7]. Under a high-frequency time-varying excitation field, the SPIOs in the FFP or FFL are dynamically magnetized and induce a voltage signal [8], [9].
To recover the SPIO concentration distribution from the measured voltage signal, the system matrix (SM) [10], [11] and x-space methods [12], [13], [14] are the main reconstruction approaches. Without prior calibration, the x-space method is less time-consuming compared with the SM-based method. X-space theory is derived under an adiabatic assumption, and the Langevin function is used to model SPIO magnetization [12], [15]. Several studies have shown that the Langevin function is inaccurate in modeling the dynamic magnetization behavior of SPIOs [15], [16], [17]. This negatively affects the spatial resolution of the reconstructed image and limits the further application of x-space MPI technology in the medical field.
The reason for the inaccuracy of the Langevin function is that it ignores the relaxation effect and assumes that SPIOs instantaneously follow the MPI magnetic field [15]. However, under the sinusoidal excitation field in MPI, SPIOs experience complex rotations, and their magnetization falls behind the excitation field, leading to the magnetic relaxation effect of SPIOs [18]. The dynamic magnetization curve of the SPIOs is a hysteresis loop in MPI, rather than the single curve described by the Langevin model [18], [19], [20], [21], [22]. Accurately modeling the dynamic magnetization behavior is essential for MPI voltage analysis and the resolution of x-space reconstruction results.
Recently, researchers have shown an increased interest in introducing appropriate models for describing the dynamic magnetization of magnetic nanoparticles in MPI to improve the image resolution [15], [19], [23]. The Debye model simplifies the magnetization behavior as a relaxation decay process [24], [25]. The Debye model uses the relaxation time constant, τ , to represent all rotations in the exponential decay process [15]. However, a single relaxation time constant cannot completely describe the rotation mechanisms and dynamic magnetization [19]. Other studies [23], [26], [27] have also investigated the exact physical causes of magnetic relaxation but have not yet applied the models to MPI reconstruction for image quality improvement. Löwa et al. proposed requirements for a dynamic magnetization model, in which the magnetization should be a function of the characterizations of SPIOs, external magnetic field, and surrounding environment [28]. From an imaging perspective, the main concern is to develop an accurate magnetization model that can be utilized to improve the spatial resolution of x-space MPI reconstruction.
The Debye model simplifies the magnetization behavior as a relaxation decay process [24], [25]. The Debye model uses the relaxation time constant, τ , to represent all rotations in the exponential decay process [15]. However, a single relaxation time constant cannot completely describe the rotation mechanisms and dynamic magnetization [19]. Other studies [23], [26], [27] have also investigated the exact physical causes of magnetic relaxation but have not yet applied the models to MPI reconstruction for image quality improvement. From an image perspective, the main concern is having an accurate magnetization model that can be utilized to improve the spatial resolution of x-space reconstruction.
To establish a more accurate magnetization model for SPIOs in MPI, the Jiles-Atherton (JA) model, which is an advanced model [29], is introduced. The JA model has good performance in modeling ferromagnetic hysteresis mechanisms and is widely utilized in numerous applications, including rotating electric machines [30], power transformers [31], and magnetic field sensors [32]. Although the JA model exhibits high accuracy when describing the magnetization of solid magnetic materials, directly applying it to the SPIOs in MPI would cause errors due to inappropriate parameter definitions. Some modifications are needed so that a more accurate magnetization model can be applied to the x-space reconstruction algorithm to improve the spatial resolution of MPI technology.
In this study, a modified Jiles-Atherton (MJA) model is developed to describe the dynamic magnetization behavior of SPIOs in MPI. Three modifications are made based on the JA model and MJA model parameters are adjustable based on the magnetization measurements. The interaction between the magnetization of SPIOs and the external fields is also included in the MJA model, making the theory closer to the magnetization process of the SPIOs in MPI. Then, the MJA model is integrated into the x-space algorithm to reduce the image blur caused by the relaxation effect and improve the image resolution. Furthermore, the precision of the MJA model is validated by comparing it with actual magnetization measurement data. The performance of the MJA x-space algorithm is evaluated by conducting simulations and one-dimensional (1D) MPI experiments.
The remainder of this paper is organized as follows. Section II presents the theory and methodology of the MJA model and the MJA x-space reconstruction algorithm. In Section III, we introduce the experiments for the MJA model and the MJA x-space reconstruction algorithm validation. The results are presented and discussed in Sections IV and V, respectively. Finally, Section VI concludes the paper.

A. X-Space MPI Theory
First, we review the signal-generating process of the xspace MPI theory. In the 1D MPI system, the gradient field is −G[A/m 2 ] and the time-varying excitation field is H 0 (t) = Acos(2πf t); f and A are the excitation frequency and amplitude, respectively. The total magnetic field is H( is the position of FFP. The magnetization of SPIOs M in response to the applied magnetic field obeys the following magnetization function: which can be described by different models, such as the Langevin theory, Debye model, and other magnetization models [15], [33]. Via the reciprocity theorem [34], [35], the induced voltage is given by where ρ ( x) denotes the concentration distribution of the SPIOs, B 1 denotes the receive coil sensitivity, * represents the convolution operator,ẋ s (t) is the velocity of the FFP, and M (x) = dM dx is the point spread function (PSF) [12], [14]. To recover the concentration of the particles from the voltage signal, the xspace algorithm first normalizes the voltage signal based on the velocity of the FFP [12] IMG native = s(t) To obtain a better image quality, the native image, IMG native , can be deconvoluted via the PSF [14] as follows: where * represents the deconvolution operator, and IMG f inal is the final image with better quality.
In the existing MPI theory, the Langevin function is a widely used magnetization function, and (1) can be written as where β := μ 0 m k B T , k B is the Boltzmann constant, m[A · m 2 ] is the magnetic moment of the SPIO, T is the particle temperature, and L(·) is the Langevin function [1].
The use of the Langevin function to model the magnetization of SPIOs under static or quasi-static magnetic fields is reasonable [36]. However, in MPI, the excitation frequency is typically up to kilohertz in order, and the Langevin function is not suitable for the SPIOs. When the sinusoidal excitation field increases and decreases, the measured magnetization curves and the derivative of the magnetization curves are intrinsically different from those of the Langevin function [20], [37], [38], as shown in Fig. 1. Using the Langevin function in MPI leads to inaccuracies in the voltage signal simulation and image blurring in the x-space reconstruction.
To address this problem, the Debye model was introduced to the MPI theory, describing the magnetization behavior as a relaxation process [15]. The Debye model adds a temporal convolution term to (7) as follows: where τ is relaxation time constant and u(t) is the Heaviside function. The Debye x-space algorithm utilizes half the relaxation time for velocity compensation However, a single relaxation time constant τ cannot completely describe the SPIO rotation mechanisms and dynamic magnetization [19]. Although using half of the relaxation time constant for velocity compensation yielded satisfactory results, it does not mean this is the best compensation value. If a more accurate compensation parameter can be derived through theoretical analysis, the image quality would be further improved.

B. MJA Model
Accurate modeling of the dynamic magnetization behavior of SPIOs is a topic of interest to the MPI community. Accordingly in this study, we propose an MJA model to describe the magnetization of SPIOs more accurately.
First, we provide a brief overview of the JA model before discussing the MJA model. The JA model can be conveniently written as a single ordinary differential equation describing the magnetization of sample M under external magnetic field H [29]: where parameter M an is the anhysteretic magnetization, and the anhysteretic function can be selected quite arbitrarily (as in [39]), α is the interdomain coupling in the magnetic material, k is the average energy required to break the pinning site in the magnetic material, and c ∈ [0, 1] is the magnetization reversibility. The JA model is usually used for solid magnetic materials [29]. For SPIOs in the liquid phase, some definitions in the JA model cannot be directly applied to the SPIOs, and some parameters are unsuited.
Based on the JA model and combined with specific SPIO features, the MJA model is adopted to accurately model the magnetization behavior of SPIOs in MPI. Next, the details of the MJA model are introduced. The MJA model for MPI is derived by adopting three modifications.
Modification One: Using coupling parameter α to express the complicated interaction between the SPIOs and the excitation field, the total effective magnetic field, H e , is defined as When an external field, H, is applied to the SPIOs, they are magnetized. Then, the magnetized SPIOs react to the magnetic field and contribute to the total effective magnetic field. This complicated interaction between the SPIOs and the excitation field is simplified and represented by α. This coupling parameter is not only determined by the characteristics of the SPIOs and the excitation field, but is also affected by the surrounding environmental variables, such as temperature, viscosity, and velocity. The determination of α is presented in Section II-C. Such an effective field expression has been used by D.C. Jiles for presenting interdomain coupling of ferromagnetic hysteresis [29].
Modification Two: Applying the Langevin function as the anhysteretic function for M an , and choosing the anhysteretic function parameters from the magnetization parameters of SPIOs.
The relationship between H e and the anhysteretic magnetization, M an , can be written as where β = μ 0 m k B T is a parameter characterizing the magnetization feature of SPIO, similar to parameter β in Section II.A, and M s is the saturation magnetization of the SPIO. Although there has been an expression similar to the formula (12) in [29], the parameter β in [29] is with dimensions of magnetic field to characterize the shape of magnetization, which is substantially different from that in the MJA model.
Modification Three: The reversible parameter c is assumed as zero to simplify the model.
If there is no simplification, the magnetization M MJA consists of two parts where the reversible magnetization, M rev , is related to the domain wall bending and the irreversible magnetization, M irr , is related to the domain wall displacement [40].
In MPI, the SPIOs align along the external magnetic field, leading to domain wall displacement. Thus, M irr is dominant, and M rev and c can be assumed as zero to simplify the model. A similar simplification has been used to model a specimen of Fe-C hysteresis magnetization curve [29].
Under the assumption that c = 0, the MJA magnetization is equal to M irr , i.e., Based on the above modifications, the final MJA model is obtained as follows: where k[A/m] is related to the coercive point; the determination of k is presented in Section II-C. Equations (11) and (12) are substituted into (16) to solve the ordinary differential equation; thus, the relationship between M MJA and H can be determined as M = M MJA (H).

C. Determination of the MJA Parameters
To obtain a more accurate and robust model, the MJA model parameters can be adjusted according to the magnetic characteristics parameters of the SPIOs. In the MJA model, some nanoparticle parameters, such as the diameter and saturation magnetization, can be easily found in the SPIO datasheet. The other two important parameters, α and k, can be determined using actual magnetization curve measurement data. The hysteresis loop of SPIOs has two important points: the coercive (H = H c , M = 0) and remanence (H = 0, M = M r ) points, as shown in Fig. 1(a).
At the coercive point with δ = +1, the MJA model equation, i.e., (16), can be written as . (18) Similarly, at the remanence point with δ = −1, (16) can be written as Then, parameters k and α can be determined by combining (18) and (19). For various excitation fields and different SPIOs, k and α vary and must be recalculated using this method. Subsequently, all parameters in the MJA model can be obtained, and the magnetization of SPIOs in MPI can be simulated.

D. X-Space MPI Theory Based on the MJA Model
Based on the MJA model, magnetization function (1) in the MPI theory can be modified with When the applied magnetic strength is H hy , the derivative of the magnetization dM dH reaches its maximum, as shown in Fig. 1(b). This value can be calculated through the MJA model, describing the extent to which the real magnetization of the SPIOs lags behind the external magnetic field.
The induced voltage in (3) can be written as The dynamic magnetization behavior causes an out-ofsynchronization between the FFP movement and the induced voltage. The voltage signal, s MJA (t), corresponds to velocitẏ x s (t − t d ) rather thanẋ s (t). The time delay t d can be calculated as follows: Integrating the MJA model into x-space reconstruction, the MJA x-space algorithm equation can be rewritten as where t d is used to compensate for the mismatch between the FFP velocity and the voltage signal. This compensation can reduce the adverse impact of dynamic magnetization on the image, resulting in a resolution improvement. For further image quality improvement, the PSF of the MJA model is used in the deconvolution process where * represents the deconvolution operator, PSF MJA is close to the actual PSF in MPI. The width of the dynamic range of the magnetization characteristic can be alternatively reported using the full-width at half-maximum (FWHM) of the derivative of magnetization. Thus, the FWHM of the MJA model is given by an empirical formula (25) where H F W HM L is the FWHM of the Langevin model. To test the accuracy of this equation, an experiment could be found in the Supplementary Section C. Because the Langevin model neglects the dynamic magnetization of SPIOs during the MPI scan, the FWHM predicted by it is inaccurate. Conversely, the MJA model takes into account the real-world magnetization phenomena; thus, the FWHM calculated by the MJA model is closer to the real FWHM of the SPIOs.

III. EXPERIMENTS
There are two main parts in this section, namely, the MJA model validation experiments (Section III-B) and the MJA x-space algorithm validation experiments (Section III-C). To validate the MJA model, the magnetization curve calculated via the MJA model was compared with the measurement data obtained via magnetic particle spectrometry (MPS). The robustness of the MJA model was tested under various test conditions, such as different excitation frequencies, amplitudes, SPIO types, and biasing fields. To test the MJA x-space algorithm, 1D simulations were conducted to examine its performance in the presence of noise. Accordingly, 1D MPI phantom experiments were performed to evaluate its image resolution improvement ability.

A. Evaluation Indicators
We used different indicators for evaluating the model magnetization describing performance and reconstructed image quality.
For a quantitative analysis of the proposed model performance in terms of describing the dynamic magnetization of SPIOs, the root mean square error (RMSE) was applied to calculate the error between the model simulated data and the measured data: where M real is the magnetization measured by the MPS, N is the number of data, and M model is the magnetization simulated by the model. To quantify the reconstruction results, the FWHM and signalto-background ratio (SBR) of the images were introduced to evaluate the reconstructed image results.
The FWHM of the reconstructed image is defined as the width at which a function decays to 50% of its maximum [36]. For the MPI reconstructed image, the FWHM of the 1D profile was calculated. A lower FWHM indicates a better spatial resolution.
The SBR is defined as the mean signal divided by the mean background to quantify the MPI reconstruction quality, as derived in [41]. The true value of the SPIOs concentration distribution was selected as a criterion for distinguishing between the signal and background areas. A higher SBR indicates better image quality.
where pixel s and pixel b represent the pixel values of the signal and noise areas, respectively.

1) Magnetic Particle Spectrometry and Nanoparticles:
To explore the dynamic magnetization behavior of SPIOs, the MPS developed by our team was used to record the magnetization curves of the nanoparticles under various sinusoidal excitation fields. The signal of the excitation coil was generated on a digital-to-analog converter, amplified, and band-pass filtered.
The induced voltage signal in the receiver coil was first canceled out by the direct coupling component from the excitation coil and then amplified and digitized by using an analog-to-digital converter. To make the exact signal timing, the voltage signal was calibrated and synchronized by the same clock with the excitation signal in the converter. Because the voltage signal is the derivative of the magnetization, the calibrated voltage signal is integrated over time to obtain the magnetization. The excitation amplitude was adjustable in the range of 0-10 mTμ −1 0 . Three excitation frequencies, namely, 1, 10, and 25 kHz, were selected, and the sampling frequency was 1 MHz. A Helmholtz coil is added perpendicular to the excitation coil, and a direct current is supplied to provide a biasing field in the range of 0-6 mTμ −1 0 . In this study, four commercial magnetic particles were selected for testing: Perimag (Micromod GmbH, Germany) coated with dextran (surface: NH2), Synomag (Micromod GmbH, Germany) coated with dextran (surface: NH2), VivoTrax (Magnetic Insight, USA) coated with carboxydextran (surface: plain), and Mag3300 (NanoEast, China) coated with oleic acid (surface: DSPE-PEG-NH2). The iron concentrations in Perimag, Synomag, VivoTrax, and Mag3300 are 5, 6, 5.5, and 1 mg/ml, respectively. Perimag consists of cross-linked dextran iron oxide composite particles with a hydrodynamic diameter of 130 nm. Synomag is available with a hydrodynamic diameter of 70 nm. VivoTrax has an average core size of 4.2 nm and a hydrodynamic size of 62 nm. Mag3300 is a single-core sample with a particle size of 20 nm. The samples were diluted to various concentrations to investigate their magnetization behavior.
2) Experimental Setup: The dynamic magnetization of the SPIOs can be affected by several factors, such as the excitation frequency, amplitude, SPIO type, and biasing field. To validate the MJA model and test its robustness under various conditions, five groups of magnetization experiments were conducted. In each group, a volume of 60 μl sample was first measured by using the MPS instrument, and the measured data were regarded as the true value. Subsequently, the MJA, Langevin, and Debye models were used to simulate the magnetization under these test conditions. The MJA model parameters were obtained by using the method described in Section II-C. The Langevin model parameters were calculated by using the method described in [12]. The relaxation time, τ , of the Debye model, was calculated by using the fitting algorithm presented in [15]. The MJA and Debye model parameters are listed in Supplementary material Tables SI-SV. The RMSE and RE values of the three models were calculated.
Group 1: Frequency test. Excitation frequencies of 1, 10, and 25 kHz were chosen for this test. The amplitude of the excitation field was 5 mTμ −1 0 and the sample was Perimag with 0.2 mg/ml iron concentration.
Group 2: Amplitude test. Excitation amplitudes of 6, 8, and 10 mTμ −1 0 were selected for this test. The excitation frequency was 10 kHz and the sample was Perimag with 0.1 mg/ml iron concentration.
The hysteresis parameter, H hy , was also calculated under these test conditions. The influence of different factors on the magnetic hysteresis effect was also explored.
Experiments in which the sinusoidal excitation field was substituted with a trapezoidal excitation field and a triangular excitation field were also conducted. However, the MJA model was originally designed for sinusoidal excitation and its performance for other excitations was not perfect. The results of these experiments are presented in the Supplementary Material.

1) 1D MPI Simulations:
One-dimensional MPI simulation experiments were performed using MATLAB. The gradient field strength was 2.5 T/m, and the excitation field was 25 mT/μ 0 at 25 kHz. The sampling frequency was set as 1 MHz. The MJA model was used to generate the voltage signal. To test the robustness of the MJA x-space algorithm under the presence of noise with a signal-to-noise ratio of 5, 10, and 20 dB was added to the simulated voltage signal. Then, three reconstruction algorithms, namely, x-space [12], Debye x-space [15], and MJA x-space, were used to reconstruct the images. The parameters of the simulations were β = 22.5 × 10 −3 A/m, k = 5.68 × 10 −3 A/m, α = 5 × 10 −4 , τ = 2 μs. The value of τ was chosen through a parameter test, which was described in the Supplementary Section D. The results were evaluated by using the SBR and FWHM.

2) 1D MPI Scanner and Phantom Experiments:
The 1D MPI scanner has a 2.6 T /m main field gradient and a 10 mTμ −1 0 excitation field at 25 kHz, as shown in Fig. 2. The sampling frequency was 1 MHz. Two cuboid-resolution phantoms with an inner thickness of 2 mm and a scaled load platform were created by using a 3D-printer. A Perimag solution with an iron concentration of 2.5 mg/ml was injected into the two cuboid phantoms. The positions of the two phantoms can be adjusted  arbitrarily, and the inner side distances between them were set to 6, 3, and 2 mm. For these experiments, three algorithms, namely, the xspace [12], Debye x-space [15], and MJA x-space, were used to reconstruct the image, and the results are evaluated using the SBR and FWHM. The time delay, t d = 2.03 μs, in the MJA x-space algorithm and the relaxation time, τ = 3.05 μs, in the Debye x-space algorithm had been calculated in the MPS experiments, as listed in Supplementary Table SI.

A. MJA Model Validation Results
From the MJA model validation results for the five experimental groups, the accuracy and wide applicability of the MJA model were demonstrated. Further, how each factor influences the dynamic magnetization lag behind the external magnetic field was also investigated based on these experiments.  The RMSE of the three models is shown in Fig. 4. For the five experimental groups, the average RMSE of the Langevin model (0.316 ± 0.09) was the largest among the three models. The average RMSE of the Debye and MJA models was 0.131 ± 0.036 and 0.055 ± 0.014, respectively. When modeling the dynamic magnetization of the SPIOs, the average RMSE of the MJA model was 83% and 58% lower than that of the Langevin and Debye models, respectively.
The robustness of the MJA model was also demonstrated in these five experimental groups. When the SPIOs are excited by sinusoidal magnetic fields with different frequencies or amplitudes, the MJA model can accurately describe the SPIOs' magnetization response. When the SPIOs experience different biasing fields in addition to the dynamic excitation field, the MJA model still describes the magnetization behavior with small errors. This is important for the MJA model applicability in MPI as the spatial encoding depends on the biasing field that forms the FFP, and the SPIOs slightly off the FFP (having a small biasing field of 0 to 5 mT μ −1 0 ) still have some contribution to the MPI signal. For different kinds of commercial SPIOs at any concentration, regardless of being single-(such as Mag3300) or multi-core (such as Perimag), the performance of the MJA model is satisfactory. All these results demonstrate that the MJA model has the ability to accurately describe the dynamic magnetization of the SPIOs under various conditions in MPI.

2) Factors Influencing the Dynamic Magnetization Lagging Extent:
Based on the MJA model, we also calculated the parameter H hy under different test conditions to explore the factors influencing the dynamic magnetization lagging extent, as shown in Fig. 5. Fig. 5(a) shows that the magnetization lagging extent is independent of the particle concentration. This is a critical point that provides a basis for introducing the hysteresis parameter, H hy , and its corresponding time delay, t d , into the MPI reconstruction algorithm, which does not generate additional errors in the concentration reconstruction results.
The magnetization lagging extent is sensitive to the excitation frequency. Fig. 5(b) shows the magnetization curves for Perimag calculated by employing the MJA model at 1, 10, and 25 kHz. Both the maximum magnetization and coercive points are easily affected by the frequency. Regardless of the type of nanoparticles, the higher the frequency, the larger the parameter, H hy , as shown in Fig. 5(c). This guides the selection of a suitable MPI excitation frequency. To reduce the magnetization lagging in MPI, a lower excitation frequency would be a better choice. However, to use the magnetization lagging feature of different SPIOs to perform multicolor or multifunction imaging, a high excitation frequency may be required.
The magnetization lagging extent increases as the excitation amplitude increases, as shown in Fig. 5(d). Because a high-amplitude excitation field is required to generate the SPIO nonlinear response in MPI, the magnetization lagging extent is too large to be ignored. Therefore, in the MPI reconstruction process, image blur caused by the magnetization lagging must be eliminated.
Finally, for different types of SPIOs, the frequency, amplitude, and iron concentration have similar effects on the magnetization lagging extent. Among the tested SPIOs, Synomag exhibited the largest magnetization lag.

1) 1-D Simulation Results:
The MJA x-space algorithm can reconstruct the SPIO concentration distribution at the presence of noise. Fig. 6(a) and (b) show the normalized SPIO concentration distribution and FFP trajectory, respectively. The pure voltage signal and the voltage with a 10 dB signal-tonoise ratio are shown in Fig. 6(c). The reconstructed 1D images obtained by the MJA x-space algorithm are shown in Fig. 6(d) when a 20, 10, and 5 dB Gaussian noise is added to the pure voltage signal. The MJA x-space algorithm still performs at a high level in the presence of Gaussian noise above 5 dB.
Then, two steps of the MJA x-space reconstruction algorithm, namely, using t d for velocity compensation and deconvolution with PSF, were separately applied to elucidate their individual contribution to the imaging performance. The native image before deconvolution and the final image after deconvolution are shown in Fig. 7. Without any compensation, the native x-space image has two enveloped peaks on each side of the true value of the SPIO distribution, which negatively affects the image resolution. With compensation, the velocity mismatch was better compensated by the MJA x-space algorithm using parameter t d than the Debye x-space algorithm using relaxation time τ , as shown in Fig. 7(a). The FWHM of the MJA x-space result is approximately 40 % less than that of the Debye x-space, and half that of the x-space reconstructed result. This means the first step can significantly improve the image resolution. The deconvolution operation also improved the image resolution in  the three algorithms. In each algorithm, the PSF was obtained based on each forward magnetization model. The calculation process of PSF referred to the derivation in [12]. Considering the scanning directions, the final PSF used for deconvolution is the sum of positive and negative scanning PSFs, both in the MJA and the Debye x-space algorithms. In the MJA x-space algorithm, the FWHM decreased by 33 % via deconvolution. Therefore, these two MJA x-space reconstruction algorithm steps both contribute to improving the reconstructed image resolution.
The MJA x-space algorithm achieve better performance than the other two algorithms in SBR and FWHM. The SBR and FWHM of the reconstructed results are summarized in Table I. Compared with the x-space method, the image resolution was improved 69% by the MJA x-space. Compared with the Debye x-space algorithm, the MJA x-space improved 10% image resolution.
2) 1D MPI Phantom Experiment Results: Compared with the x-space and Debye x-space methods, the MJA x-space method exhibited the best performance in reconstructing the 1D MPI image and improving the image resolution. For better visualization, the 1D images were expanded from 1 × 100 to 10 × 100 pixels. Fig. 8 presents the normalized reconstructed images and the 1-D profile through the center of the images. The SBR and FWHM evaluation indicators are listed in Table II. FWHM value larger than the center-to-center distance of the phantoms means that the reconstruction process failed and the phantoms cannot be clearly distinguished in the image.
When the phantoms are placed at intervals of 6 and 3 mm, the image resolution is improved by the MJA x-space algorithm. For a 2 mm interval, the MJA x-space method can clearly distinguish the two phantoms, whereas the other methods reconstruct an incorrect phantom with a larger width. The FWHM of the x-space and Debye x-space reconstructed images is significantly larger than the center-to-center distance of the phantoms, 3 mm; therefore, these two methods failed to reconstruct a correct image. The MJA x-space algorithm can compensate for image blur better than the Debye x-space algorithm. The image resolution of the MJA x-space algorithm is, on average, 64% and 48% better than the x-space and Debye x-space algorithms, respectively.
These results indicate that the MJA x-space method can reduce the image blur caused by the dynamic magnetization behavior and improve the resolution of MPI reconstructed images.

V. DISCUSSION
In current x-space MPI, because of the existence of the relaxation effect, the Langevin model is inaccurate for describing the dynamic magnetization of SPIOs and leads to a degraded spatial resolution [15], [17]. To address this problem, in this study, we proposed an accurate and robust MJA model to describe the dynamic magnetization of SPIOs and improve the x-space image resolution in MPI. Through an ordinary differential equation, the relaxation effect of the dynamic magnetization curve of the SPIOs in MPI can be well expressed in the MJA model. Three modifications contribute to the MJA model being accurate and robust in describing various SPIO magnetizations under different excitation conditions. By integrating the MJA model into the x-space theory, the mismatch between the velocity and the MPI voltage signal is well compensated for, and the reconstructed image resolution is improved.
The good performance of the MJA model in terms of accuracy and robustness is due to the following reasons. First, the ordinary differential equation coming from the JA model is an excellent function to describe the hysteresis loop curve. Thus, it is suitable for describing dynamic magnetization with the presence of the relaxation effect. Then, with the modifications, the interaction between the SPIOs and the external field is considered in the MJA model, making the model closer to the real SPIO magnetization process in MPI. When the excitation field changes or the surrounding SPIO environment changes, the interaction is affected, and these alterations will be reflected in the MJA model. Thereafter, the MJA model parameters are determined in combination with the characteristics of the SPIOs and excitation conditions, making the model adjustable for describing different test conditions. The magnetization in the MJA model is regarded as a function of the SPIOs, external magnetic field, and surrounding environment, satisfying the requirements for a dynamic magnetization model proposed in [28]. Based on these advantages, the MJA model has the ability to model the dynamic magnetization behavior with high accuracy and robustness.
Precisely because the MJA model is accurate, it is effective to integrate the MJA model into the x-space algorithm and use the time delay, t d , to reduce the image blur and improve its resolution. Even though there is possible to obtain a rough value of time delay directly from the MPS data, the accuracy is worse than that of the theoretical value calculated by the MJA model. There could be a relationship between t d in the MJA model and the relaxation time, τ , in the Debye model, which will be investigated in future work.
In this work, we compared the MJA x-space algorithm with the x-space and Debye x-space methods, which improve the MPI reconstruction by changing the magnetization model. Regarding state-of-the-art x-space reconstruction algorithms, such as those in [42] and [43], we did not consider the reconstruction quality improvement based on aspects other than the magnetization model. The MJA model, as a more accurate magnetization model, could potentially be combined with these advanced x-space algorithms to further improve image quality.
Despite these advantages, the MJA model also has some limitations. Although the MJA model performs well in describing the magnetization under sinusoidal excitation fields, for the other excitation waveforms, its performance suffers somewhat. The magnetization curves calculated by the MJA model under trapezoidal and triangular excitation fields are shown in Figure S1 (see the Supplementary Material). A potential improvement could be achieved by applying more useful points in the previously measured magnetization curve to calculate the MJA model parameters. Furthermore, compared with the Debye model with only one adjustable parameter, τ , the MJA model is more complex and has three adjustable parameters, k, α, and β, to improve its accuracy and robustness. This complexity issue could be addressed by preparing a parameter dictionary and trade-off between model complexity and accuracy. Finally, the current theoretical analysis and experimental results show that the model accuracy or image quality was not affected with the assumption of c=0 for model simplification. Thus, it is acceptable to assume c as zero. Whether using the full JA model without this assumption will benefit requires further study.
The MJA model also has the potential to be applied in simulating the system matrix. Acquiring the calibration-based system matrix is time-consuming procedure. An accurate magnetization model of the SPIOs can be used to simulate the system matrix to alleviate the time cost issue in system matrix acquisition. Several relevant research works have studied the simulation of the system matrix [44], [45], [46]. However, the accuracy of the simulated system matrix still needs for further improvement.
Additionally, the influence of the surrounding environment of SPIOs on the magnetic hysteresis effect will be explored in future work. According to the current MJA model analysis, the surrounding environment can affect the coupling parameter and further change the time delay; however, this specific relationship will be investigated with more experiments. The MJA model has the potential to be used in multiparameter (such as viscosity, temperature, and velocity) imaging [25].

VI. CONCLUSION
In this study, the MJA model is proposed to describe the dynamic magnetization behavior of SPIOs with accuracy and robustness. As a connection between the real-world dynamic magnetization and the reconstructed image, the MJA model can be applied in the x-space algorithm to improve the spatial resolution of MPI technology. Considering the relaxation effect, and the interaction between the SPIOs and the magnetic field, the MJA model becomes more accurate and robust than the Langevin and Debye models. Compared with the x-space and Debye x-space methods, the MJA x-space algorithm exhibits better performance in terms of improving the spatial resolution of the MPI reconstruction image. The MJA model also has the potential to be extended to velocity and temperature imaging in MPI.