Calibration-Free Relaxation-Based Multi-Color Magnetic Particle Imaging

TL;DR

This paper introduces a calibration-free method for multi-color magnetic particle imaging that estimates nanoparticle relaxation times directly from the signal, enabling distinction of different nanoparticles without prior calibration.

Contribution

The paper presents a novel relaxation time estimation technique based on mirror symmetry, eliminating the need for extensive calibration in multi-color MPI.

Findings

Successfully distinguishes nanoparticles without calibration.

Validated method through simulations and experiments.

Effective at multiple frequencies (550 Hz and 9.7 kHz).

Abstract

Magnetic Particle Imaging (MPI) is a novel imaging modality with important applications such as angiography, stem cell tracking, and cancer imaging. Recently, there have been efforts to increase the functionality of MPI via multi-color imaging methods that can distinguish the responses of different nanoparticles, or nanoparticles in different environmental conditions. The proposed techniques typically rely on extensive calibrations that capture the differences in the harmonic responses of the nanoparticles. In this work, we propose a method to directly estimate the relaxation time constant of the nanoparticles from the MPI signal, which is then used to generate a multi-color relaxation map. The technique is based on the underlying mirror symmetry of the adiabatic MPI signal when the same region is scanned back and forth. We validate the proposed method via extensive simulations, and via experiments on our in-house Magnetic Particle Spectrometer (MPS) setup at 550 Hz and our in-house MPI scanner at 9.7 kHz. Our results show that nanoparticles can be successfully distinguished with the proposed technique, without any calibration or prior knowledge about the nanoparticles.

Equations

$$s_{adiab}(t)=B_{1}m\rho (x)\ast \dot{\mathcal{L}}[kGx]{|}_{x=x_s(t)}kG{\dot{x}}_s(t)$$

$$s(t)=s_{adiab}(t)\ast r(t)$$

$$r(t)=\frac{1}{\tau }{e}^{-\frac{t}{\tau }}u(t)$$

$$s_{pos,adiab}(t)=-s_{neg,adiab}(-t)=s_{half}(t)$$

$$s_{pos}(t)=s_{pos,adiab}(t)\ast r(t)=s_{half}(t)\ast r(t)$$

$$s_{neg}(t)=s_{neg,adiab}(t)\ast r(t)=-s_{half}(-t)\ast r(t)$$

$$\mathcal{F}\{r(t)\}=R(f)=\frac{1}{1+i2\pi f\tau }$$

$$\mathcal{F}\{s_{pos}(t)\}=S_{pos}(f)=S_{half}(f).R(f)$$

$$\mathcal{F}\{s_{neg}(t)\}=S_{neg}(f)=-S_{half}^{\ast }(f).R(f)$$

$$\hat{\tau }(f)=\frac{S_{pos}^{\ast }(f)+S_{neg}(f)}{i2\pi f(S_{pos}^{\ast }(f)-S_{neg}(f))}$$

$$\hat{\tau }=\frac{{\int }_0^{f_{max}}|S_{pos}(f)|\tau (f)df}{{\int }_0^{f_{max}}|S_{pos}(f)|df}$$

$$({\hat{t}}_{edge},\hat{\tau })=\underset{(t_{edge},\tau )}{\mathrm{argmin}}{\int }_0^{T/2}({\hat{s}}_{pos}(t)-(-{\hat{s}}_{neg}(-t){)}^{2}dt$$

$$s_{pos}(t)=c_1{s}_{half}(t)\ast r_1(t)+c_2{s}_{half}(t)\ast r_2(t)$$

$$s_{neg}(t)=c_1(-s_{half}(-t)\ast r_1(t))+c_2(-s_{half}(-t)\ast r_2(t))$$

$$R_i(f)=\frac{1}{1+i2\pi f{\tau }_i},i=1,2$$

$$S_{pos}(f)=S_{half}(f).(c_1{R}_1(f)+c_2{R}_2(f))$$

$$S_{neg}(f)=-S_{half}^{\ast }(f).(c_1{R}_1(f)+c_2{R}_2(f))$$

$$\hat{\tau }(f)$$

$$=\frac{2c_1{\tau }_1(1+4{\pi }^2{f}^2{\tau }_2^2)+2c_2{\tau }_2(1+4{\pi }^2{f}^2{\tau }_1^2)}{2c_1(1+4{\pi }^2{f}^2{\tau }_2^2)+2c_2(1+4{\pi }^2{f}^2{\tau }_1^2)}$$

$$\hat{\tau }(f)\approx \frac{c_1{\tau }_1+c_2{\tau }_2}{c_1+c_2}$$

Full text

I Introduction

Magnetic Particle Imaging (MPI) is a new and rapidly developing imaging modality, that shows great potential in terms of resolution, contrast, and sensitivity [1, 2, 3, 4, 5]. This imaging technique exploits the nonlinear magnetization response of superparamagnetic iron oxide (SPIO) nanoparticles which provides MPI with high-contrast imaging capabilities, without any signal from the human tissue [1, 6, 7]. MPI uses combination of three different magnetic fields to obtain an image: a static selection field with a strong gradient is used to create a field-free-point (FFP). Then a sinusoidal drive field is applied to move the FFP and scan the field-of-view (FOV). However, human safety limits restrict the size of the FOV covered by the drive field [8, 9, 10, 11, 12]. Hence, low-frequency focus fields are used to shift the FFP to regions that cannot be covered with the small FOV due to drive field alone [13, 14, 15].

There are two main image reconstruction methods in MPI. System Function Reconstruction (SFR) requires calibration of the nanoparticle response and the imaging system, which is achieved by a time-consuming calibration scan that records the frequency response of the signal induced by the SPIOs located at all voxel locations in the FOV [16, 17, 18, 19]. The alternative image reconstruction technique, x-space reconstruction, does not require any calibration as it reconstructs the image by gridding the acquired signal to the instantaneous position of the FFP [20, 21]. The resulting images, in return, show a version of the SPIO distribution blurred by the point spread function (PSF) of the imaging system. Both of these image reconstruction techniques, with their own advantages and disadvantages, map the spatial distribution of the SPIOs. While MPI’s quantitative imaging of SPIO distributions can be an important tool by itself for applications such as cardiovascular imaging and stem cell tracking, separation of the signals acquired from different nanoparticles would increase the functionality of MPI even further.

Rahmer et al. have recently demonstrated that different SPIO types can be distinguished via a multi-color MPI technique [22], and numerous applications have already been shown to benefit from this novel approach. Recently, there has been significant progress in using multi-color MPI for catheter tracking during cardiovascular interventions [23]. In such applications, one SPIO type is injected into the blood stream for vessel visualization, while the catheter is coated with a different type of SPIO for tracking purposes. A further advancement of this technique involved simultaneous tracking and steering of the catheter tip via using the the magnetic fields in MPI [24, 25]. More recently, the temperature mapping capability of multi-color MPI has also been demonstrated [26].

Binding state of the SPIOs has been spectroscopically shown to change relaxation behavior as well [27]. In this regard, relaxation mapping can be used to probe cell or protein binding of SPIOs. Moreover, there have been efforts on drug delivery via nanocarriers [28], where multi-color MPI can offer additional tracking possibilities. Likewise, multi-color MPI can also be used to identify the characteristics of an environment such as viscosity and temperature, which are shown to affect nanoparticle relaxation behavior significantly [29, 30, 31]. These conditions can be an important tool for probing the changes in tissue environments, e.g., in the case of hyperthermia treatments. All of the aforementioned applications create an immense room for further experimental work and research.

To date, multi-color MPI has been realized via both SFR and x-space approaches. In the SFR approach, SPIOs were differentiated based on the differences in their harmonic responses, which were obtained using an extensive calibration procedure performed separately for each SPIO type [22]. For the x-space approach, Hensley et al. demonstrated that x-space reconstruction is capable of multi-color reconstruction via differentiating relaxation behaviors of different SPIO types, if multiple measurements at different drive field amplitudes are utilized [32].

In this work, we present a novel, calibration-free multi-color MPI technique for x-space MPI. This technique can generate a relaxation time constant map of SPIOs from a single scan at a single drive field amplitude. The proposed technique takes advantage of the back and forth scanning of a FOV to estimate the relaxation time directly from the MPI signal, without any a priori knowledge about the SPIOs. Here, we demonstrate with extensive simulation results, MPS experiments, and imaging experiments that the proposed method successfully distinguishes multiple SPIO types, without requiring any calibrations.

II Theory

II-A MPI Signal

In x-space MPI, the ideal signal (also known as the adiabatic signal) is defined via the Langevin response of the nanoparticles to an applied drive field [20]:

[TABLE]

Here, $B_{1}~{}[T/A]$ is receive coil sensitivity, $m~{}[Am^{2}]$ is the magnetic moment of nanoparticles, $\rho(x)~{}[particles/m^{3}]$ is nanoparticle density, $k~{}[m/A]$ is nanoparticle property, and $G~{}[T/m/\mu_{0}]$ is the gradient strength of the selection field. In practice, the MPI signal lags and gets wider due to nanoparticle relaxation effects, i.e., the delay of SPIOs in aligning with oscillating magnetic fields (see Fig. 1a). In x-space MPI, the relaxation effect is modelled as a temporal convolution of the ideal MPI signal with an exponential kernel [33]:

[TABLE]

Here, r(t) denotes the relaxation kernel and $\tau$ is the relaxation time constant. The resulting MPI signal is called the non-adiabatic signal. Via extensive experimental studies, this simple yet powerful model has been shown to accurately match the relaxation effect for a wide range of drive field frequencies and amplitudes [33, 34].

II-B Proposed Method: Direct Relaxation Time Constant Estimation

In MPI, a sinusoidal drive field is applied to excite the SPIOs, while additional focus fields are applied to control the global positioning of the FFP. For a fixed focus field amplitude, the resulting FFP trajectory scans a partial field-of-view (pFOV) back and forth around a central position. For simplicity, the forward motion of the FFP is called the positive scan, while the backward motion is called the negative scan. For that trajectory, one would expect the signals acquired during these two scans to be mirror symmetric, as displayed in Fig. 1a-b. However, the nanoparticle relaxation causes an asymmetric blurring that breaks the mirror symmetry (see Fig. 1c).

We have previously introduced a technique to directly estimate the relaxation time constant from the MPI signal, using the underlying mirror symmetry of the adiabatic signal [35, 36, 31]. In this technique, we formulated the effects of relaxation in Fourier domain to directly estimate the relaxation time constant. Here, we define $s_{pos}(t)$ as the signal acquired during positive scan, and $s_{neg}(t)$ as the signal acquired during negative scan, both centered with respect to time. For the adiabatic MPI signal, $s_{pos,adiab}(t)$ and $s_{neg,adiab}(t)$ are mirror symmetric, i.e.,

[TABLE]

where, $s_{half}(t)$ denotes half a period of the adiabatic MPI signal. Using the relaxation formulation in Eq. 2, the non-adiabatic positive and negative signals can be expressed as:

[TABLE]

Here, $s_{half}(t)$ and r(t) are the unknowns of the equations while $s_{pos}(t)$ and $s_{neg}(t)$ are measured waveforms. In Fourier domain, these equations can be expressed as:

[TABLE]

where $\mathcal{F}$ is the Fourier transform operator, and the time reversal and conjugate symmetry properties of Fourier Transform are used to express $S_{neg}(f)$. Using Eqs. 7-9 , $\tau$ can be calculated directly as follows:

[TABLE]

Ideally, performing this calculation at a single frequency in Fourier domain suffices. However, the accuracy of the estimation is different at each frequency component due to relative strengths of signal vs. noise. To increase the robustness of the estimation against noise, a weighted average of ${\tau}(f)$ is calculated with respect to the magnitude spectrum:

[TABLE]

Here, $f_{max}$ is an upper threshold for the range of frequencies used, such that $0<f_{max}\ll\frac{F_{s}}{2}$, where $F_{s}$ is the sampling frequency. Typically, including frequencies up to $6^{th}$ or $7^{th}$ harmonic of the fundamental frequency suffices, as the signal falls of rapidly with increasing frequency.

II-C Verification with 1D MPI Simulations

Figure 2 demonstrates the proposed relaxation time constant estimation method on simulated MPI signal. Here, a triangular nanoparticle density is assumed and a central pFOV is scanned to obtain the MPI signal. Then, relaxation effect is simulated via the exponential model in Eq. 2, with $\tau=7~{}{\mu}s$. Figure 2b shows the asymmetric blurring due to relaxation, for which the estimated time constant yielded $\hat{\tau}=6.96~{}{\mu}s$. The signal was then deconvolved using the estimated relaxation kernel, recovering the underlying mirror symmetry (see Fig. 2c). The same procedure was repeated on signal after direct feedthrough filtering [6, 37], to ensure that the filtering of the fundamental harmonic (a necessary step in MPI) does not hinder the performance of the proposed technique. The filtered signal is shown in Fig. 2d with the corresponding deconvolved signal in Fig. 2e. Here, the estimated time constant was $\hat{\tau}=6.99~{}{\mu}s$. In the absence of noise in both cases, the proposed method estimated $\tau$ with less than 0.6% error, which is primarily caused by the digitization in simulating the relaxation effect via convolution.

III Materials and Methods

III-A Experimental Relaxation Time Constant Estimation

Correct timing of $s_{pos}(t)$ and $s_{neg}(t)$ is crucial for the proposed relaxation time estimation method. In simulations, this procedure is trivial since the timing of the signal vs. the FFP trajectory (i.e., the time point, $t_{edge}$, when the positive scan ends and the negative scan starts) is known. In practice, however, system delays introduce a time shift between the signal vs. the FFP trajectory. In such cases, incorrect $t_{edge}$ values introduce extra phase terms in Eq. 10, causing incorrect estimation of $\tau$. One method for measuring $t_{edge}$ is performing a calibration scan with an SPIO that does not exhibit relaxation, such that the zero-crossing of the signal acquired during calibration scan can be assigned as $t_{edge}$. Here, we propose a procedure that estimates the correct $t_{edge}$ and $\hat{\tau}$ simultaneously, eliminating the need for a calibration scan.

We have argued previously that deconvolution of the MPI signal with the correct relaxation kernel should restore mirror symmetry. Here, we exploit this phenomenon to simultaneously find the correct $t_{edge}$ and $\hat{\tau}$, as outlined in Fig. 3. First, we choose a $t_{edge}$ value, directly estimate $\tau$ using Eq. 11, and deconvolve the MPI signal with the estimated relaxation kernel. We then repeat this step for a range of $t_{edge}$ values limited to [0, T/2) region of the MPI signal, where T is the period of the FFP trajectory. Correct ($t_{edge},\tau$) pair should restore the mirror symmetry, minimizing the mean squared error (MSE) between positive and mirrored negative signals after the deconvolution, i.e.,

[TABLE]

where $\hat{s}(t)$ denotes the signal after deconvolution with the estimated relaxation kernel. Here, instead of computing MSE over the entire half period, more weights can be assigned to central time points that typically have higher SNR [31].

III-B Proposed Algorithm: Calibration-Free Multi-Color MPI

For the proposed multi-color MPI technique, we directly estimate the relaxation time constant for each pFOV, and map the estimated time constant to the corresponding spatial location. While this mapping sounds simple at first, extensive simulations revealed two special cases where the additional steps are needed for accurate mapping of time constants:

Flat Nanoparticle Distribution: If the nanoparticle distribution is flat in a pFOV, most of the signal is contained in the lost 1st harmonic. Hence, direct feedthrough filtering removes almost all signal, making flat regions appear like regions devoid of nanoparticles. In the absence of signal, proposed estimation method produces noise-like results. 2. 2.

Inhomogeneous Mixtures of Different Nanoparticle: In the regions where two different nanoparticle types mix homogeneously, the estimation yields a weighted average of two relaxation time constants (see Appendix). However, in the case of an inhomogeneous mixture (e.g., a transition region from one nanoparticle type to the other), the adiabatic signal from each nanoparticle not only has a different shape, but is also convolved with a different relaxation kernel. Therefore, the set of equations derived in Eqs. 8-9 become underdetermined, producing noise-like estimations.

Here, we first detect and isolate these special cases, reconstruct the rest of the relaxation map, and finally restore the missing portions of the map. The following steps summarize the proposed multi-color MPI algorithm, outlined in Fig. 4, that can handle these special cases:

Phase 0 - Estimation of $\tau$ for Each pFOV: Using the signals from each pFOV, algorithm directly estimates the corresponding relaxation time. An x-space MPI image is also reconstructed (Fig. 4b), and the signal RMS values for pFOVs are calculated (Fig. 4c).

Phase 1 - Histogram Correction: Using the histogram of estimated $\tau$ values, a mild thresholding is applied to remove unlikely estimations The histogram correction consists of 2 steps (see Fig. 4d-e):

• •

Upper limit on estimations: Extensive work in nanoparticle relaxometry showed that relaxation time constants are much smaller than the period of the drive field [34], typically less than 10% of the period [31]. Here, we ignore estimations that are larger than one-fourth of the period.

• •

Signal RMS threshold: Low signal leads to inaccurate $\tau$ estimations. Here, estimations from low RMS regions (e.g., less than 10% of the maximum RMS value) are initially ignored.

Phase 2 - Cluster Detection: At the end of Phase 1, unlikely estimations are mostly removed. Accordingly, the remaining estimations must have the relaxation times of different SPIO types. Here, a k-means clustering algorithm is run on the remaining estimations [38]. Since the number of nanoparticle types is unknown, $k$ parameter is started from a large value and iteratively reduced until convergence is achieved. Reduction method is as follows (see Fig. 4f-g):

• •

Merger of clusters: Cluster centers that are closer than a $\tau_{res}$ value are merged. Note that $\tau_{res}$ must be sufficiently small to ensure the separation of different clusters.

• •

Omission of small clusters: Clusters that contain less than $p\%$ of the overall estimation number are considered irrelevant and ignored.

• •

Convergence: If the cluster properties (number of clusters, cluster centers, etc.) remain the same for two consecutive iterations, convergence is achieved.

Phase 3 - Special Case Detection: After Phase 2, clusters are detected, which means probable SPIO types and their relaxation time estimations are identified. Outlier $\tau$ values that are not in the vicinity of the cluster centers could indicate either one of the two special cases explained above. These cases have high pixel intensity in the MPI image despite having low signal RMS. Here, we search for the nearest spatial locations that have $\tau$ values belonging in a cluster, on either side of the problem pFOV. If incorrect estimations are caused by inhomogeneous mixtures, we expect different clusters on either sides of the problem pFOVs. On the other hand, flat nanoparticle distributions have a single cluster type around the problem pFOVs (see Fig. 4h).

Phase 4 - Recovery of $\tau$ Map: At the end of Phase 3, the regions of incorrect estimations are determined. These regions are recovered as follows (see Fig. 4i):

• •

For Flat Nanoparticle Distributions, the missing regions are recovered as a constant $\tau$ value.

• •

For Inhomogeneous Mixtures, the average signal RMS of the clusters are used to estimate the relative concentrations of the corresponding nanoparticle types. Next, a transition region from one relaxation time to the other is modeled as a sigmoid function, with steepness determined via relative concentrations.

Finally, estimations from low MPI pixel intensity regions are set to zero to remove the background in the $\tau$ map.

III-C 1D and 3D Simulations

All simulations were carried out using a custom MPI toolbox developed in MATLAB (Mathworks, Natick, MA). To obtain realistic results, fundamental harmonic was filtered out [6], relaxation time constant values were selected according to previous studies [34, 31], and smoothly-varying Hamming-window-shaped SPIO distributions were used. The regular MPI images were reconstructed using the x-space reconstruction technique [20, 21, 6].

For 1D simulations, a selection field with 3 $[T/m/\mu_{0}]$ gradient was utilized, with 10 mT-peak drive field at 10 kHz and 80% overlap between consecutive pFOVs. The nanoparticle diameter was selected as 21 nm, a conservative value given the recent developments in tailored SPIO production [39]. The relaxation times ranged between 1-5 $\mu s$, as recently reported in experimental studies performed around 10 kHz [34, 31]. To match experimental conditions, the simulated MPI signal was sampled at 2 MSPS sampling frequency, and additive white Gaussian noise corresponding to a peak signal SNR of 100 was added to the MPI signal.

For 3D simulations, a selection field with (-7, 3.5, 3.5) $[T/m/\mu_{0}]$ gradient in (x, y, z) directions was utilized, based on the work in [6]. A 15 mT-peak drive field at 10 kHz and a 90% overlap along z-direction was utilized, with all other parameters kept the same as in 1D simulations.

III-D MPS Experiments

The initial experiments of proposed relaxation time estimation method were conducted on our in-house Magnetic Particle Spectrometer (MPS, also known as MPI Relaxometer). The drive coil of the relaxometer produces 0.97 mT/A magnetic field with 95% homogeneity in a 7 cm long region. The receive coil of the relaxometer was designed as a three-section gradiometer pick-up coil [40].

Experiments were conducted at 550 Hz and 15 mT-peak drive field. The nanoparticle signal was first amplified with a low-noise voltage preamplifier (Stanford Research Systems SR560), then digitized at 2 MSPS via a data acquisition card (National Instruments, NI USB-6363). 16 consecutive acquisitions were averaged to increase the SNR. A background measurement was subtracted from the averaged signal to minimize the effect of direct feedthrough. A frequency domain filter was applied by selecting the higher harmonics of the drive field and setting the rest of the frequency components to zero in Fourier domain. Next, a high-order zero-phase digital low pass filter (LPF) was applied in time domain to restrict the range of harmonics used. The cut-off frequency of the LPF was determined by comparing noise power with signal power at the harmonic frequencies, while avoiding the 250 kHz self-resonance frequency of the receive coil. A more detailed explanation of the complete setup can be found in [31].

III-E Imaging Experiments

Our in-house MPI Scanner (see Fig. 5) has two disc-shaped permanent magnets with 2-cm thickness and 7-cm diameter, placed at 8-cm separation in x-direction. The resulting configuration creates (-4.8, 2.4, 2.4) $[T/m/\mu_{0}]$ gradient in (x, y, z) directions, which yields approximately 4-mm resolution in z-direction (i.e., down the imaging bore) for Nanomag-MIP nanoparticles. The drive coil has 3 layers of Litz wire with 80 turns, resulting in 1.5 mT/A magnetic field with 95% homogeneity in a 4.5 cm-long region down its bore. The drive field and selection field specifications were validated using a Hall Effect Gaussmeter (LakeShore 475 DSP Gaussmeter). The receive coil was designed as a three-section gradiometer, with a single layer of Litz wire with 34 turns in the main section and 17.5 turns in the side sections. The self-resonance of the receive coil was measured at around 280 kHz. This configuration allowed 1x1x10 cm3 FOV. The drive coil and receive coil were placed inside a cylindrical copper shield with 1-cm thickness that can be used for imaging at drive field frequencies as low as 1 kHz.

Figure 5a-c display the workflow of the complete MPI scanner setup, together with front and side views. Overall imaging system was controlled via a custom toolbox implemented in MATLAB (Mathworks, Natick, MA). The drive coil was impedance matched to AE Techron 7224 power amplifier using a capacitive network at 9.7 kHz. The imaging phantoms were mechanically moved inside the scanner via Motor-Driven Velmex BiSlide (Model: MN10-0100-E01-21) in 3 dimensions.

The drive field was at 15 mT-peak and 9.7 kHz, which resulted in a 12.5 mm pFOV length in z direction. The amplitude of the drive field was calibrated immediately before each experiment via a Rogowski current probe (LFR 06/6/300, PEM Ltd.). Partial FOVs were acquired with 85% overlaps. For 1D imaging experiments, an 8-cm FOV along z-direction was covered with 16.8 sec active scan time. For 2D imaging experiments, a recti-linear trajectory was used to cover a 0.8x6.8 cm2 FOV in x-z plane using 9 lines along z-direction, with 134 sec active scan time. The temperature inside the scanner bore was controlled throughout the experiment to prevent heating of the nanoparticles.

To boost the SNR, 16 consecutive acquisitions were averaged for each pFOV. The background measurements acquired before/after each line were subtracted from the nanoparticle signal to remove direct feedthrough. All remaining signal acquisition/processing steps were the same as in the MPS experiments. The regular MPI images were reconstructed using x-space reconstruction [20, 21, 6], followed by the proposed multi-color MPI technique.

III-F Phantom Preparation

Imaging phantoms with different types of SPIOs were prepared, using Vivotrax ferucarbotran nanoparticles (Magnetic Insight Inc., USA) with same chemical composition as Resovist, and Perimag and Nanomag-MIP nanoparticles (Micromod GmbH, Germany). These nanoparticles had original concentration levels of 5.5 mg Fe/mL, 17 mg Fe/mL, and 5 mg Fe/mL, respectively. To approximately equalize the MPI signal levels, Perimag and Nanomag-MIP SPIOs were diluted 10.2 and 3.5 times, respectively. Capillary tubes with 2 mm inner diameter were filled with 10 $\mu$L volume (after dilution) of SPIOs each. These tubes were then placed in custom-designed 3D-printed phantom holders in two different configurations, as shown in Fig. 5d-e.

IV Results

IV-A 1D and 3D Simulation Results

Figure 6 displays the results of the proposed mapping algorithm for the 1D case, where each row corresponds to a different scenario: two, three, and four different types of SPIOs, respectively. The left column shows the regular MPI images for SPIO with relaxation times ranging between 1-5 $\mu s$. The overall MPI image is a summation of these images (not displayed for clarity). Each simulation includes a challenging case of inhomogeneous mixtures of different nanoparticle types. The right column shows the reconstructed relaxation maps (blue solid lines) and ideal relaxation maps (orange dashed lines). While utilizing only two different SPIO types is more likely in medical applications (e.g., catheter tip vs. blood pool in the case of cardiovascular interventions [25]), these results aim to show the full potential of the proposed method.

As seen in each row in Fig. 6, proposed mapping algorithm reconstructed the relaxation maps accurately, despite the presence of noise. The inhomogeneous mixture regions were detected successfully and the proposed sigmoid-transition model recovered those regions accurately. Furthermore, the proposed algorithm also eliminated background areas via a thresholding and edge detection step. The slight deviations from the ideal relaxation maps were due to noise, as well as the crosstalk of signals from different SPIOs due to proximity. For the results in Fig. 6b, the mean estimation error is well below 3%, whereas the mean error reaches 7% for the results in Fig. 6f where the SPIO distributions are closer to each other. Here, the 3.5-$\mu$s SPIO was spatially disconnected from the other SPIOs, which is successfully depicted in the reconstructed relaxation map.

For the multidimensional case, the proposed technique is extended in a line-by-line basis, where each line in z-direction is reconstructed individually, then combined to form the multidimensional multi-color MPI image. In Fig. 7a, each SPIO distribution (labeled from 1 to 9) represents a homogeneous mixture of two different SPIOs with 2.9 $\mu$s and 1.1 $\mu$s relaxation time constants at different mixture ratios. Accordingly, the distribution labeled as 1 has 100% of 2.9-$\mu$s SPIOs, whereas the distribution labeled as 9 has 100% of 1.1-$\mu$s SPIOs. Figure 7b-c show the resulting regular MPI image and the color overlay of relaxation map, respectively. In the regular MPI image, the regions corresponding to larger $\tau$ resulted in slightly lower pixel intensities due to relaxation-induced signal loss. However, this difference cannot be used to distinguish the SPIOs, as lower iron concentration could also result in reduced pixel intensity. As seen in Fig. 7c, the proposed multi-color MPI method is capable of distinguishing a variety of SPIO types (9 in this case), each with a different time constant. Furthermore, the estimated time constants reflect the relative concentrations of the constituent SPIO types with high level of linearity ($R^{2}=0.9921$), as shown in Fig. 7d.

IV-B MPS Experiment Results

Figure 8c provides proof-of-concept validation for the simulations given in Fig. 7. In the MPS experiments, nanomag-MIP particles were mixed with water/glycerol at two different ratios, to obtain two different relaxation times due to differences in viscosity levels [31]. As given in Fig. 8c, relaxation time constants for these two cases were estimated as 64.83 $\mu$s and 46.22 $\mu$s.

Next, a homogeneous mixture was mimicked by acquiring signal from both samples simultaneously. From theory (see Appendix), the resulting relaxation time constant should be the average of the relaxation time constants of the constituents weighted by their concentrations. Here, since the two samples had the same concentration of SPIOs, the expected relaxation time of the homogeneous mixture is 55.52 $\mu$s. The estimated relaxation time of 54.21 $\mu$s closely mathces this expectation.

IV-C Imaging Experiment Results

Figure 9 displays the 1D imaging results of the proposed technique. In this experiment, a phantom with 3 tubes, each containing a different type of SPIO was imaged. To enable color display, reconstructed 1D regular MPI image and 1D relaxation map were replicated and stacked vertically in a pseudo-2D image format. Figure 9a shows the reconstructed regular MPI image. Although the SPIOs were placed in identical tubes, Vivotrax covers a wider region (7.3 mm) in the MPI image, denoting a lower resolution capability. On the other hand, Perimag and Nanomag-MIP SPIOS display comparable widths (5 mm vs. 5.2 mm).

While there are differences between MPI responses of these SPIOs, it is not possible to distinguish the SPIO types based on the regular MPI image in Fig. 9a. Figure 9b shows the relaxation map reconstructed with the proposed algorithm. Here, different SPIO types can be clearly distinguished after color assignment. On average, Nanomag-MIP, Perimag, and Vivotrax SPIOs yielded 2.7 $\mu$s, 3.5 $\mu$s and 4.5 $\mu$s relaxation time constants, respectively. Finally, Fig. 9c is the color overlay image, containing the spatial and quantitative information from both the regular MPI image and the relaxation map.

Figure 10 displays 2D experimental demonstration of the proposed method. Here, Nanomag-MIP and Vivotrax SPIOs were imaged with a 12 mm separation. Again, while the two SPIOs were placed in identical tubes, Vivotrax displays a significantly wider PSF in both x- and z-directions in the regular MPI image in Fig. 10a. In both the relaxation map and the color overlay image, Nanomag-MIP and Vivotrax SPIOs can be distinguished clearly. The average time constants were measured as 2.5 $\mu$s and 4.3 $\mu$s, respectively. Note that these results closely match the relaxation times measured in the 1D experiments for Nanomag-MIP and Vivotrax SPIOs. Here, the colormap (shown in horizontal bar) was chosen to be identical to that in the 1D experiments, to enable a direct comparison of results. Accordingly, 1D and 2D results show a strong agreement, indicating a capability for quantitative imaging using the proposed multi-color MPI technique.

In Fig. 10b, at the peripheries of the SPIO distributions along the x-direction, the estimated relaxation times deviate with respect to those in more central regions. While this is an undesired effect of low signal, color overlay image naturally eliminates these low-pixel-intensity regions to provide a clean multi-color MPI image, as seen in Fig. 10c.

V Discussion

Multi-color MPI is a very promising extension of MPI due to its potential use in cardiovascular interventions, temperature mapping, and viscosity mapping. There has been significant progress in multi-color MPI, both via SFR and x-space reconstruction techniques. While these methods were shown to work exceptionally well, they either require extensive calibration for every SPIO types used [22] or multiple scans of the imaging FOV at different drive field amplitudes [32]. In this work, we have proposed and experimentally demonstrated a calibration-free approach for multi-color MPI, that works successfully with a single scan at a single drive field amplitude. Here, calibration-free refers to SPIO identification through their relaxation times, which are estimated directly from the MPI signal without any prior information. Note that while relaxation times can be measured in a calibration-free fashion, for certain applications such as viscosity or temperature mapping, one must obtain prior fiducial measurements to create a dictionary of relaxation times corresponding to specific environmental conditions. Alternatively, these measurements could be obtained directly during imaging via placing the fiducials next to the object of interest.

Distinguishable SPIO behavior is a fundamental requirement for any multi-color MPI technique. Rahmer et al. have suggested that, SPIOs with similar size distributions or similar hysteresis behavior would be difficult to separate [22]. A similar requirement applies in our case, where SPIOs need to have distinct relaxation time constants. For the experimental results shown in this work, all three nanoparticles involved were multi-core SPIOs, with Perimag and Nanomag-MIP having similar chemical compositions and hydrodynamic diameters. Despite these similarities, the proposed multi-color MPI technique successfully distinguished these nanoparticles. To further improve the performance of this technique, two different approaches can be undertaken: SPIOs that are specifically tailored to have different relaxation behaviors would facilitate relaxation mapping. In addition, imaging parameters can be optimized to maximize the contrast between relaxation times. Our previous work on relaxation-based viscosity mapping showed that relatively low drive field amplitudes and frequencies (e.g., 1 kHz) are more suitable for viscosity mapping, provided that the SNR is sufficiently high [31]. On the other hand, for distinguishing different SPIO types, it may be more advantageous to operate at higher frequencies and drive field amplitudes. A previous work has compared the relaxation times of multi-core Resovist nanoparticles with single-domain UW33 nanoparticles, which had significantly different relaxation times. Interestingly, the ratio of the relaxation times remained almost constant (around 2-2.5) across a wide range of frequencies and drive field amplitudes [34]. In such a case, operating at higher fields/frequencies would be desirable for improving the SNR of the MPI signal, although the human safety limits of the drive field should also be taken into account [9, 8, 11].

In this work, we implemented a conventional FFP trajectory, where a piecewise-constant focus field was mimicked with mechanical movement of the robot arm in the x- and z-directions, while a sinusoidal drive field repetitively scanned a chosen pFOV in z-direction. This resulted in 134 sec active imaging time for the 2D case. The scan times can be significantly shortened by utilizing more efficient FFP trajectories. For example, a high-slew-rate linearly-ramping focus field can be utilized in place of a piecewise-constant focus field. Our preliminary simulation results for such cases yield accurate relaxation time estimations for slew rates as high as 10 T/s (results not shown). With such a trajectory, 134 sec imaging time would reduce down to approximately 150 ms. Experimental demonstration of rapid imaging with the proposed technique remains a future work.

VI Conclusion

In this work, we have proposed a calibration-free multi-color MPI technique and demonstrated the results via both MPS and imaging experiments. The proposed technique relies on distinguishing nanoparticles based on their relaxation response, where the relaxation time is directly estimated from the MPI signal via restoring the underlying mirror symmetry. A step-by-step algorithm is proposed to generate an accurate relaxation map (or multi-color MPI image) from the time constant estimations for each pFOV. The imaging results show that different nanoparticle responses can be successfully distinguished, without any calibration or prior information regarding their responses. The proposed calibration-free multi-color MPI technique is a promising method for future functional imaging applications of MPI, such as catheter tracking, viscosity mapping, temperature mapping, and stem cell tracking.

[Relaxation Time Constant of a Homogeneous Mixture] As stated in Section III-B, when two different SPIOs are mixed homogeneously, $\hat{\tau}$ is the weighted average of the relaxation time constants of the two SPIOs. Here, we assume that the mixture is completely homogeneous and the two SPIOs have identical adiabatic signals, $s_{half}(t)$ (e.g., bound/unbound SPIOs with identical nanoparticle core sizes). Therefore, Eqs. 5 and 6 can be modified as follows:

[TABLE]

Here, $c_{1}$ and $c_{2}$ are the concentrations of the two different SPIOs, and $r_{1}(t)$ and $r_{2}(t)$ are the corresponding relaxation kernels with $\tau_{1}\neq\tau_{2}$. Fourier domain representation of Eqs. 13 and 14 can be written as follows:

[TABLE]

Here, there are three unknowns, $S_{half},R_{1}(f)$ and $R_{2}(f)$, and two equations, Eq. 16 and 17. Hence, the equation set above is underdetermined. Nevertheless, if the original estimation method given in Eq. 10 is applied, the following expression can be obtained using Eqs. 16-17.:

[TABLE]

Equation VI is not strictly linear with respect to concentrations of the SPIOs in the mixture. However, for small frequencies (e.g., for $f<4f_{0}$ where $f_{0}$ is the drive field frequency), the terms $(1+4\pi^{2}f^{2}\tau_{i}^{2})$ can be approximated as 1. In that case, Eq. VI simplifies to the following form:

[TABLE]

where $\hat{\tau}(f)$ is equal to the weighted average of $\tau_{1}$ and $\tau_{2}$.

Note that, if the adiabatic responses of the two SPIOs are different (e.g., SPIOs with different core sizes), the estimated $\hat{\tau}(f)$ will be an average of $\tau_{1}$ and $\tau_{2}$, weighted by the overall signal levels from the two SPIOs.

Acknowledgment

A preliminary version of this work was presented at the $6^{th}$ International Workshop on Magnetic Particle Imaging (IWMPI 2016). This work was supported by the Scientific and Technological Research Council of Turkey (TUBITAK 114E167), by the European Commission through FP7 Marie Curie Career Integration Grant (PCIG13-GA-2013-618834), by TUBA-GEBIP 2015 program of the Turkish Academy of Sciences, and by BAGEP 2016 Award of the Science Academy.