2D Quantitative Imaging of Magnetic Nanoparticles by an AC Biosusceptometry Based Scanning Approach and Inverse Problem

The use of magnetic nanoparticles (MNPs) in biomedical applications requires the quantitative knowledge of their quantitative distribution within the body. AC Biosusceptometry (ACB) is a biomagnetic technique recently employed to detect MNPs in vivo by measuring the MNPs response when exposed to an alternate magnetic field. The ACB technique presents some interesting characteristics: non-invasiveness, low operational cost, high portability, and no need for magnetic shielding. ACB conventional methods until now provided only qualitative information about the MNPs’ mapping in small animals. We present a theoretical model and experimentally demonstrate the feasibility of ACB reconstructing 2D quantitative images of MNPs’ distributions. We employed an ACB single-channel scanning approach, measuring at 361 sensor positions, to reconstruct MNPs’ spatial distributions. For this, we established a discrete forward problem and solved the ACB system’s inverse problem. Thus, we were able to determine the positions and quantities of MNPs in a field of view of 5×5×1 cm3 with good precision and accuracy. The results show the ACB system’s capabilities to reconstruct the quantitative spatial distribution of MNPs with a spatial resolution better than 1 cm, and a sensitivity of 1.17 mg of MNPs fixed in gypsum. These results show the system’s potential for biomedical application of MNPs in several studies, for example, electrochemical-functionalized MNPs for cancer cell targeting, quantitative sensing, and possibly in vivo imaging.


Introduction
Magnetic nanoparticles (MNPs) present several desirable characteristics for diagnostic and therapeutic applications, where the MNPs can be employed for tumor detection, heat generation, or drug delivery [1,2]. Therefore, techniques capable of detecting and quantifying MNPs' distribution within specific body regions, such as in tumors and organs, are essential to optimize and evaluate biomedical procedures using MNPs [3,4].
Magnetorelaxometry (MRX), Magnetic Particle Imaging (MPI), and Magnetic Susceptibility Imaging (MSI) are promising techniques proven to detect, quantify, and reconstruct spatial distributions based on the MNPs' magnetic properties. Each of these techniques is based on a specific property of the MNPs. The MRX relies on the time-dependent relaxation not straightforward. The calculations are commonly performed by solving elliptic integrals, which can involve a significant computational cost and extensive time. More recently, Hanson and Hirshman (2002) introduced a simple numerical method of the Biot-Savart law, enabling faster calculations of magnetic fields produced by filamentary segments [25]. Here, we adapted the Batuscheck and Williamson formulation, and employed the Hanson and Hirshman method and Faraday's law, to establish a forward model for the ACB system, which enabled both the discretization of the problem (i.e., establishing voxels) and linking the MNPs' mass in each voxel to the ACB signal acquired. By minimizing the least square differences between model and measurements, the most probable locations and quantities of MNPs causing the detected signals are found as the solution of the inverse problem [23].
Due to the problem characteristics (i.e., low information and a high number of possible solutions), the problem is ill-conditioned. Several methods have been proposed to improve the conditioning of inverse problems. Liebl and co-workers (2014) experimentally demonstrated that the use of multiple inhomogeneous magnetizing fields can improve the inverse problem conditioning in MRX, thus presenting better results than conventional homogeneous fields [26]. Moreover, Baumgarten and Haueisen (2010) showed that including temporal information in the MRX model also improved the inverse problem conditioning [27]. In general, the addition of information in the forward model (e.g., increasing the number of equations in the problem) improves the system matrix conditioning and stabilizes the solution of the inverse problem [26], allowing more accurate quantifications and a higher spatial resolution to be obtained.
Previously, Jaufenthaler et al. (2020) solved an ill-conditioned inverse problem to perform quantitative 2D MRX imaging of MNPs by employing optically pumped magnetometers (OPMs) [4]. The group showed that, by using OPMs and solving an inverse problem, it is possible to develop a portable MNP imaging setup. In this study, we aimed to develop an ACB setup for 2D quantitative imaging of MNP distributions, which is envisioned for ex vivo quantitative imaging of MNP distributions in organs and tissues. For instance, Wiekhorst et al. (2012) dissected pig lungs into smaller measurable pieces to quantify the MNP distribution in the whole organ using MRX [3]. Employing the presented setup, it is possible to obtain this quantitative information in a non-destructive and automatized manner. At the same time, the knowledge gained in terms of technology and methodology will serve as a starting point to develop advanced imaging setups towards 3D and in vivo quantitative MNP imaging by ACB.
Here, we aimed to develop for the first time an ACB scanning approach, and a discrete ACB forward model and its inverse solution, to quantitatively image 2D distributions of MNPs. This methodology is able to determine the MNPs' mass and position with good precision and accuracy, and much better spatial resolution (at least 1 cm) than in the previous ACB studies. The scanning approach allows obtaining a high sensitivity (reaching a minimum of 1.17 mg of MNPs) due to short sensor-to-source and excitation coilto-source distances (0.5 cm) over the entire FOV and an adjustable spatial resolution that depends on the defined scan grid. Compared to a multi-channel ACB imaging approach, the required acquisition time is much higher. However, this drawback is not critical for ex vivo samples. Gypsum phantoms with different concentrations of MNPs were considered to study the ability of the ACB system to perform quantitative imaging of MNPs. In addition, we compared both experimental and calculated system matrices for solving the inverse problem and reconstructing the MNP distributions in the phantoms. Additionally, the results presented here are also applicable to other ACB imaging modalities and will be of great value for the development of new imaging setups.

ACB Single-Channel System
The ACB single-channel system consists of two pairs of excitation and pick-up coils, in which the detection coils are connected in a first-order gradiometric arrangement, as shown in Figure 1. The system works as a double magnetic flux transformer with an air core. It can be divided into two sub-systems: the measurement and the reference system. Both subsystems are composed of one excitation and one pick-up coil, with a coaxial arrangement where the inner coil is the pick-up coil. In practice, the measurement system is positioned near the sample, and is responsible for exciting and detecting the MNPs. The reference system is placed far from the sample because of its long baseline (i.e., it is not sensitive and does not contribute to the magnetization of MNPs), detecting the environmental noise and the exciting magnetic field contributions, which will be subtracted by the gradiometric arrangement. Theoretically, no signal arises in the pick-up coils when no MNPs are near the measurement system. In practice, the system is balanced to reach a common mode rejection of at least 1.2 × 10 −4 . When an MNP sample is near the sensor, there is an imbalance in the total system's magnetic flux, and an electric signal is generated in the pick-up coils.
reference system is placed far from the sample because of its long baseline (i.e., it is not sensitive and does not contribute to the magnetization of MNPs), detecting the environmental noise and the exciting magnetic field contributions, which will be subtracted by the gradiometric arrangement. Theoretically, no signal arises in the pick-up coils when no MNPs are near the measurement system. In practice, the system is balanced to reach a common mode rejection of at least 1.2 × 10 −4 . When an MNP sample is near the sensor, there is an imbalance in the total system's magnetic flux, and an electric signal is generated in the pick-up coils.
We used a lock-in and an audio amplifier to apply a current in the excitation coils, which will generate the excitation field for the magnetizing of the MNPs. As a result of the pick-up coils' gradiometric configuration, the lock-in amplifier offers signal pre-processing and reduces environmental noise, removing the need for magnetic shielding.
We built the single-channel sensor ( Figure 1) with pick-up coils with an internal diameter of 10.14 mm, an external diameter of 11.84 mm, and 450 turns (AWG 36), and excitation coils with an inner diameter of 13.50 mm, an external diameter of 22.3 mm, and 150 turns (AWG 24). The coils were wound in a multi-layer's solenoid configuration with a width of 10 mm, and a baseline of 15 cm between the pair of the measurement and reference coils.

ACB System Forward Model and Reconstruction
According to Bastuscheck and Williamson (1995), the magnetic flux generated in the pick-up coils by an extensive magnetized source can be described as: We used a lock-in and an audio amplifier to apply a current in the excitation coils, which will generate the excitation field for the magnetizing of the MNPs. As a result of the pick-up coils' gradiometric configuration, the lock-in amplifier offers signal pre-processing and reduces environmental noise, removing the need for magnetic shielding.
We built the single-channel sensor ( Figure 1) with pick-up coils with an internal diameter of 10.14 mm, an external diameter of 11.84 mm, and 450 turns (AWG 36), and excitation coils with an inner diameter of 13.50 mm, an external diameter of 22.3 mm, and 150 turns (AWG 24). The coils were wound in a multi-layer's solenoid configuration with a width of 10 mm, and a baseline of 15 cm between the pair of the measurement and reference coils.

ACB System Forward Model and Reconstruction
According to Bastuscheck and Williamson (1995), the magnetic flux generated in the pick-up coils by an extensive magnetized source can be described as: where φ s is the magnetic flux generated in the pick-up coil, µ 0 is the magnetic permeability in a vacuum, Ir is the pick-up coil reciprocal current, χ the material's magnetic susceptibility, B a the excitation field, and B r the reciprocal field [24]. This model was conveniently built using the reciprocity principle, where we can interchange detector and source positions to compute the induced signal. Using the reciprocity principle allows application of the numerically efficient expression of the Biot-Savart law by Hanson and Hirshman (2002) for both excitation and sensing [25]. By Hanson and Hirshman's model, the magnetic flux density, at the center of k-th voxel, produced by a coil, is given by where p 1,i,k and p 2,i,k are the vectors that respectively connect the start point and the final point of i-th segment of the coil's wire to center of the k-th voxel. Using Equation (2), each turn of the coil can be approximated by 36 rectilinear segments to maintain the calculations errors below 1% [26]. For quantitative 2D imaging, a discrete forward model for ACB is necessary. Thus, we discretize the problem, and consequently, the FOV is composed of K voxels. Hence, the voltage induced in the pick-up coil, generated by an MNPs mass (X MNP,k ) in the k-th voxel, can be calculated by Faraday's law and has a final form as: where V n,k is the voltage induced in the pick-up coil, H r,n,k and H e,n,k are the reciprocal and the magnetizing field, respectively, at the k-th voxel with the n-th position sensor. Note that in our ACB setup, the MNPs are magnetized by an alternate magnetic field with a constant frequency. Thus, the frequency-dependent mass susceptibility χ(ω) of the MNPs becomes constant in our model. Equation (3) can be extended to N sensor positions and K voxels. Using matrix notation and separating the geometric parameter and MNPs' properties from the MNPs' mass within each voxel, we obtain the following expression: where V is a vector containing the induced voltages in each sensor, L is the sensitivity matrix of the system, which comprises both geometric parameters and MNPs' properties, and contains all the sensitivity coefficients that link the mass of MNPs in the k-th voxel to a measurement signal in the sensor at the n-th position, and X MNP is the vector that represents the MNP mass in each voxel. It is important to point out that the sensitivity matrix L can be determined a priori by calculations or measurements, as described later. Additionally, considering a system with N sensors (positions) and K voxels will result in a matrix L ∈ R N×K . The inverse problem in ACB is often ill-posed; thus, multiple MNP distributions may explain the measurement data. To find the most probable MNP distribution in our imaging problem, we applied linear optimization techniques, such as minimum-norm estimation, employing the Moore-Penrose pseudo-inverse (L + = L T L −1 L T ) calculated by Truncated Singular Value Decomposition (TSVD), as has been shown to be feasible for similar imaging problems [26,28]:X whereX MNP is a vector containing the estimated MNPs' mass in each voxel. This work determined a truncate limit of 5% for the experimental matrix and 1% for the calculated matrix.

MNP Phantom Development
We employed custom-made MNPs, with a manganese ferrite core (MnFe2O4), synthesized by the coprecipitation method as described previously [29], with a mean diameter of 15 ± 5 nm measured using Transmission Electron Microscopy, and saturation of magnetization of 49.4 emu/g, measured by Vibrating Sample Magnetometry. Fe and Mn content were around 74.4% and 25.6%, respectively, as shown previously. Further details of the MNPs' synthesis and characterization can be found in [18,30].
Different quantities of MNPs were immobilized in gypsum for the experimental measurements. Thus, the MNPs' Brownian relaxation was suppressed and Néel relaxation was the dominant relaxation mechanism (so that the effective relaxation time equal approximately to apply the Néel relaxation time τ e f f ≈ τ N ). By drying the MNP gypsum mixture in defined containers, cubic MNP phantoms with a volume of 1 cm 3 were constructed. Five cubes with 0.43, 1.61, 6.77, 9.68, 12.38, and 12.88 mg (A, B, C, D, E1, and E2, respectively) of MNPs were constructed for the 1 cube distribution sensibility test and the 2 cubes distribution spatial resolution test. Moreover, were made 10 similar cubes with an MNP mass of approximately 6 mg to perform several cubes' distributions in the geometry of the letters: B, I, O, M, A and G.

Experimental and Calculated Sensitivity Matrices Acquisition
To perform quantitative imaging of MNP distributions, we employed a scanning approach, performed by a 2D CNC (Computed Numeric Control) stage, which moved the sensor in a 9 × 9 cm 2 grid with a step of 0.5 cm. The scanning approach represents N = 361 sensor positions in the forward model, used to acquire the MNP signals in 5 × 5 × 1 cm 3 FOV, with 1 cm 3 voxels and centralized with a sensor grid. Figure 2 shows an example of the ACB single-channel scanning of an MNP distribution in a 2D voxel grid.

MNP Phantom Development
We employed custom-made MNPs, with a manganese ferrite core (MnFe2O4), synthesized by the coprecipitation method as described previously [29], with a mean diameter of 15 ± 5 nm measured using Transmission Electron Microscopy, and saturation of magnetization of 49.4 emu/g, measured by Vibrating Sample Magnetometry. Fe and Mn content were around 74.4% and 25.6%, respectively, as shown previously. Further details of the MNPs' synthesis and characterization can be found in [18,30].
Different quantities of MNPs were immobilized in gypsum for the experimental measurements. Thus, the MNPs' Brownian relaxation was suppressed and Néel relaxation was the dominant relaxation mechanism (so that the effective relaxation time equal approximately to apply the Néel relaxation time ≈ ). By drying the MNP gypsum mixture in defined containers, cubic MNP phantoms with a volume of 1 cm 3 were constructed. Five cubes with 0.43, 1.61, 6.77, 9.68, 12.38, and 12.88 mg (A, B, C, D, E1, and E2, respectively) of MNPs were constructed for the 1 cube distribution sensibility test and the 2 cubes distribution spatial resolution test. Moreover, were made 10 similar cubes with an MNP mass of approximately 6 mg to perform several cubes' distributions in the geometry of the letters: B, I, O, M, A and G.

Experimental Sensitivity Matrix ( )
To perform quantitative imaging of MNP distributions, we employed a scanning approach, performed by a 2D CNC (Computed Numeric Control) stage, which moved the sensor in a 9 × 9 cm 2 grid with a step of 0.5 cm. The scanning approach represents = 361 sensor positions in the forward model, used to acquire the MNP signals in 5 × 5 × 1 cm 3 FOV, with 1 cm 3 voxels and centralized with a sensor grid. Figure 2 shows an example of the ACB single-channel scanning of an MNP distribution in a 2D voxel grid. During experimental signal acquisition, a lock-in amplifier and an audio amplifier applied a root-mean-square (RMS) current of 500 mA oscillating at 10 kHz in the excitation coils. The same lock-in measured the voltage ( ) induced in the detection coils for each of the N sensor positions.
We built a custom-made controlling system to perform all measurements in the Lab-View ® and Mach3Mill ® environment. Briefly, after being converted to the voltage RMS value by the lock-in amplifier (with a time constant of 1 ms), the signal was recorded at a sampling frequency of 20 Hz using an A/D acquisition board (National Instruments-NI-DAQPad 6015). Simultaneously, the CNC stage moved the sensor, stopping 0.5 s for the During experimental signal acquisition, a lock-in amplifier and an audio amplifier applied a root-mean-square (RMS) current of 500 mA oscillating at 10 kHz in the excitation coils. The same lock-in measured the voltage (V n ) induced in the detection coils for each of the N sensor positions.
We built a custom-made controlling system to perform all measurements in the LabView ® and Mach3Mill ® environment. Briefly, after being converted to the voltage RMS value by the lock-in amplifier (with a time constant of 1 ms), the signal was recorded at a sampling frequency of 20 Hz using an A/D acquisition board (National Instruments-NIDAQPad 6015). Simultaneously, the CNC stage moved the sensor, stopping 0.5 s for the voltage acquisition in each position. Furthermore, the CNC stage generated a clock signal with a width of 0.5 s (acquisition period in each position) collected in the same A/D board.
The experimental sensitivity matrix (L exp ) was built performing 25 scans, in which the sensor was moved in each of the 361 positions (Figure 2a). Each scan was performed with the cube D1 positioned in one of the K = 25 voxels. Therefore, we obtained a sensitivity matrix L exp ∈ R 361×25 . Furthermore, we built an experimental sensitivity map that summarized each sensor position's contributions in each voxel. All measurements were made with a distance of 5 mm between the sensor and cube surfaces.
Similarly, the calculated sensitivity matrix L cal was constructed through the computational simulation of the ACB system, as shown in Figure 2b.
Because the magnetic flux calculation does not have a simple analytic solution, especially in non-symmetric situations, the numerical model introduced in Equation (2) was applied to determine H r,n,k and H e,n,k . Two considerations were made for using this method: (i) each turn of the coils was polygonally composed of 36 segments, and (ii) the magnetizing field at one point was the sum of fields from all segments.
For the calculations, the real values of the current amplitude (0.5 A) and frequency (10 kHz) were considered. Regarding L exp , L cal was built with the virtual sensor positioned in each of the 361 positions. A single voxel with MNP was considered to occupy a different virtual voxel in each scanning, thus resulting in L cal ∈ R 361×25 . All the calculations were conducted considering a 5 mm distance between the voxel and sensor surfaces. Once again, we built a sensitivity map using L cal by summing the signals of all sensor positions for each cube position in a voxel.

Experiments and Quantifications
Sensitivity analysis measurements were assessed scanning five MNPs' cubes separately (cubes A, B, C, D, and E2), where each cube was positioned in the central voxel of the FOV. We further scanned two similar cubes (E1 and E2), placing the cubes with distances between the cube's surfaces of 0, 1, and 3 cm along the FOV's central line. Using the reconstructed images' correlation and quantification absolute relative error, as detailed below, we determined a threshold distance (spatial resolution in centimetric voxel) above which the system can quantitatively separate the two cubes.
Furthermore, similar cubes of 6 mg were placed and scanned in the FOV to compose B, I, O, M, A, G letters. These data were used to check the method performance in reconstructing three or more cubes at the same time, and the influence of the number of scanning points in the reconstructed image, and to qualitatively compare the voltage maps (no inverse problem solution) with the quantitative images (inverse problem solution).
To assess the influence of the number of scanning points (i.e., the number of coils in the forward problem), the letters' distributions were reconstructed and evaluated with a horizontal and vertical grid step of 0.5, 1.0, and 1.5 cm, to obtain 361, 100, and 49 grid points, respectively. All of these matrices presented the same voxel, scanning grid, and FOV geometry and sizes.
For each data set, we applied a minimum norm reconstruction based on L exp and L cal , respectively. For imaging quality evaluation, all reconstructed distributionsX MNP were correlated with the nominal distribution X nom using the Pearson Correlation Coefficient (CC) [31,32]: where CC = 1 is obtained when nominal and reconstructed distribution are equal, and CC = 0 when there is no statistical similarity between both distributions. For quantification precision assessments, the absolute relative MNPs' mass difference, X di f f , was calculated by comparing the sum of all voxel values in the reconstructed image and the nominal distribution by Equation (7): We also studied the dependence of the CC and X di f f as a function of the MNPs' mass to evaluate the impact of the signal-to-noise ratio (SNR) on the imaging quality.

Calculated Voltage Vectors
To evaluate the system performance in reconstructing smaller voxels, we calculated a L matrix with 361 scanning points in a FOV of 5 × 5 × 1 cm 3 and reduced the voxel length to 0.5 × 0.5 × 1 cm 3 , for a total of 100 voxels. The calculated signal was obtained by solving the forward problem (Equation (4)) and adding noise based on real noise spectra, as proposed by Coene et al. (2014) [33]. The calculated signals were obtained from a distribution of bar phantoms with 5 × 0.5 × 1 cm 3 and 50 mg of MNPs. Two bars were positioned with distances between the bars' surfaces of 0, 5 and 10 mm. In addition, four bars were crossed to qualitatively evaluate the spatial resolution. Furthermore, one bar was positioned in the central line of the FOV, and the line spread function (LSF) was obtained to quantitatively determine the spatial resolution by the full width at half maximum (FWHM) method [34].

Results
The experimental and calculated sensitivity maps are shown in Figure 3. The results showed that the sensitivity map built from L cal presented greater uniformity when compared to the sensitivity map built from L exp . A low correlation (CC = 0.38) was obtained by comparing both sensitivity maps, which was expected due to the noise inherent in the L exp measurements. Furthermore, the average pixel intensity of the maps is 3.13 ± 0.28 mV/mg for L exp , and 2.81 ± 0.01 mV/mg for L cal . These differences in variances and mean values can be related to the real system's (L exp ) intrinsic noise, experimental inaccuracies (CNC stage positioning and vibration), and the approximations used in the computational mode. where = 1 is obtained when nominal and reconstructed distribution are equal, and = 0 when there is no statistical similarity between both distributions. For quantification precision assessments, the absolute relative MNPs' mass difference, , was calculated by comparing the sum of all voxel values in the reconstructed image and the nominal distribution by Equation (7): We also studied the dependence of the and as a function of the MNPs' mass to evaluate the impact of the signal-to-noise ratio (SNR) on the imaging quality.

Calculated Voltage Vectors
To evaluate the system performance in reconstructing smaller voxels, we calculated a matrix with 361 scanning points in a FOV of 5 × 5 × 1 cm 3 and reduced the voxel length to 0.5 × 0.5 × 1 cm 3 , for a total of 100 voxels. The calculated signal was obtained by solving the forward problem (Equation (4)) and adding noise based on real noise spectra, as proposed by Coene et al. (2014) [33]. The calculated signals were obtained from a distribution of bar phantoms with 5 × 0.5 × 1 cm 3 and 50 mg of MNPs. Two bars were positioned with distances between the bars' surfaces of 0, 5 and 10 mm. In addition, four bars were crossed to qualitatively evaluate the spatial resolution. Furthermore, one bar was positioned in the central line of the FOV, and the line spread function (LSF) was obtained to quantitatively determine the spatial resolution by the full width at half maximum (FWHM) method [34].

Results
The experimental and calculated sensitivity maps are shown in Figure 3. The results showed that the sensitivity map built from presented greater uniformity when compared to the sensitivity map built from . A low correlation ( = 0.38) was obtained by comparing both sensitivity maps, which was expected due to the noise inherent in the measurements. Furthermore, the average pixel intensity of the maps is 3.13 ± 0.28 mV/mg for , and 2.81 ± 0.01 mV/mg for . These differences in variances and mean values can be related to the real system's ( ) intrinsic noise, experimental inaccuracies (CNC stage positioning and vibration), and the approximations used in the computational mode.  Nominal and reconstructed MNP distributions for a single MNP loaded cube located in the center of the FOV are shown in Figure 4, where the MNPs' mass of the cube was varied. For image reconstruction, the Moore-Penrose inverse pseudo-matrices L + exp and L + cal were determined and used to calculate the most probable MNPs' distributions and quantities using Equation (5). All Pearson correlation coefficients and relative MNPs' mass differences between the reconstructed and nominal distributions are shown in Table 1. The cubes of higher concentrations (B, C, D, and E2) showed high correlation values and a superior performance in L cal reconstructions. The method did not accurately reconstruct the A cube for both matrices, having high CC values, probably related to system noise. Regarding the quantification accuracy, small differences X di f f were obtained for cubes C, D, and E2, with no considerable differences between L exp and L cal . For cube B, the method showed a moderate decrease in quantification accuracy, and for cube A a severe decrease in quantification was seen with no significant difference between L exp and L cal . This can be related to the fact that L cal has lower basal noise, and with this, higher effective reconstruction quality, but reconstructions by both matrices are not sufficient to quantify mass values lower than 1.17 mg (mass value obtained from fitted curves, Figure 5, that result in CC ≥ 0.7 and X di f f ≤ 25% simultaneously) of this MNP type. Furthermore, as the scanning area is 9 × 9 cm 2 and the FOV area is 5 × 5 cm 2 , we did not find differences in the 1 cube reconstruction in the periphery region (not shown). Nominal and reconstructed MNP distributions for a single MNP loaded cube located in the center of the FOV are shown in Figure 4, where the MNPs' mass of the cube was varied. For image reconstruction, the Moore-Penrose inverse pseudo-matrices + and + were determined and used to calculate the most probable MNPs' distributions and quantities using Equation (5). All Pearson correlation coefficients and relative MNPs' mass differences between the reconstructed and nominal distributions are shown in Table 1. The cubes of higher concentrations (B, C, D, and E2) showed high correlation values and a superior performance in reconstructions. The method did not accurately reconstruct the A cube for both matrices, having high CC values, probably related to system noise. Regarding the quantification accuracy, small differences were obtained for cubes C, D, and E2, with no considerable differences between and . For cube B, the method showed a moderate decrease in quantification accuracy, and for cube A a severe decrease in quantification was seen with no significant difference between and . This can be related to the fact that has lower basal noise, and with this, higher effective reconstruction quality, but reconstructions by both matrices are not sufficient to quantify mass values lower than 1.17 mg (mass value obtained from fitted curves, Figure 5, that result in ≥ 0.7 and ≤ 25% simultaneously) of this MNP type. Furthermore, as the scanning area is 9 × 9 cm 2 and the FOV area is 5 × 5 cm 2 , we did not find differences in the 1 cube reconstruction in the periphery region (not shown).     The nominal and reconstructed MNP distributions for two MNP cubes (E1 and E2), measured simultaneously in the FOV, are presented in Figure 6. Using both and , we were able to quantitatively resolve two cubes separated by 1 and 3 cm, where each cube is visible as an individual MNPs' source in the reconstruction. The correlation coefficients for these distances were high for all cube distances, with a higher performance in reconstruction (Table 2). Moreover, the quantification showed good accuracy with no significant differences between the and methods. For the reconstruction of two cubes placed next to each other (i.e., 0 cm distance), a good correlation ( > 0.98) and good quantification accuracy ( < 5%) were obtained. It is worth pointing out that each cube occupies half of the central voxel for 0 cm of distance, summing their contributions in this voxel, and half of the next voxel in the line, contributing with half of their MNPs' mass in these voxels. These results indicate that was no difference between and results for 2 cubes' reconstruction, and using this method, the system resolution is at least 1 cm. The nominal and reconstructed MNP distributions for two MNP cubes (E1 and E2), measured simultaneously in the FOV, are presented in Figure 6. Using both L exp and L cal , we were able to quantitatively resolve two cubes separated by 1 and 3 cm, where each cube is visible as an individual MNPs' source in the reconstruction. The correlation coefficients for these distances were high for all cube distances, with a higher performance in L cal reconstruction (Table 2). Moreover, the quantification showed good accuracy with no significant differences between the L exp and L cal methods. For the reconstruction of two cubes placed next to each other (i.e., 0 cm distance), a good correlation (CC > 0.98) and good quantification accuracy (X di f f < 5%) were obtained. It is worth pointing out that each cube occupies half of the central voxel for 0 cm of distance, summing their contributions in this voxel, and half of the next voxel in the line, contributing with half of their MNPs' mass in these voxels. These results indicate that was no difference between L exp and L cal results for 2 cubes' reconstruction, and using this method, the system resolution is at least 1 cm. ; and the third row shows reconstructions using . The columns indicate the distance between both cubes, e.g., for cube E1 and E2 with the face-to-face distance of 0, 1, and 3 cm. The reconstructions of calculated voltage vectors related to the bar distributions are shown in Figure 7, and the correlation and quantification parameters are shown in Table  3. In these reconstructions, the voxels were reduced from 1 × 1 × 1 cm 3 to 0.5 × 0.5 × 1 cm 3 , for a total of 100 voxels in the same FOV of 5 × 5 × 1 cm 3 . As and showed similar results in the reconstruction of 25 voxels, these results were obtained only by the calculated method. The CC values indicate that is no difference between horizontal and vertical reconstructions. Moreover, the method was able to simultaneously solve the horizontal and vertical bars, with a slight decrease in geometry reconstruction performance. In terms of quantification, all geometries showed high precision ( < 1%) with no significant difference between each other. These results show that the scanning method has at least 5 mm of spatial resolution. Using one bar in the vertical disposition, we obtained the LSF of the image to quantitatively evaluate the resolution (image not shown). This reconstructed bar showed a of 0.95 and of 1.42%, and the FWHM was 4.59 mm. Nominal and reconstructed MNPs' distributions with two MNP loaded cubes located in the FOV at the same time. The first row is the nominal distributions; the second row shows the reconstructions using L exp ; and the third row shows reconstructions using L cal . The columns indicate the distance between both cubes, e.g., for cube E1 and E2 with the face-to-face distance of 0, 1, and 3 cm. The reconstructions of calculated voltage vectors related to the bar distributions are shown in Figure 7, and the correlation and quantification parameters are shown in Table 3. In these reconstructions, the voxels were reduced from 1 × 1 × 1 cm 3 to 0.5 × 0.5 × 1 cm 3 , for a total of 100 voxels in the same FOV of 5 × 5 × 1 cm 3 . As L exp and L cal showed similar results in the reconstruction of 25 voxels, these results were obtained only by the calculated method. The CC values indicate that is no difference between horizontal and vertical reconstructions. Moreover, the method was able to simultaneously solve the horizontal and vertical bars, with a slight decrease in geometry reconstruction performance. In terms of quantification, all geometries showed high precision (X di f f < 1%) with no significant difference between each other. These results show that the scanning method has at least 5 mm of spatial resolution. Using one bar in the vertical disposition, we obtained the LSF of the image to quantitatively evaluate the resolution (image not shown). This reconstructed bar showed a CC of 0.95 and X di f f of 1.42%, and the FWHM was 4.59 mm. 021, 21, x FOR PEER REVIEW 12 of 18 Figure 7. Reconstructed distributions of MNP loaded bars' calculated signal. Note that lower voxels were used in the same FOV of 5 × 5 × 1 cm 3 , for a total of 100 voxels of 0.5 × 0.5 × 1 cm 3 . The first row shows the vertical distributions; the second row shows the horizontal distributions; and the third row shows the cross distributions. The columns indicate the distance between the bars' surfaces. The reconstruction of the MNPs' distributions of several cubes arranged to form the letters B, I, O, M, A, and G can be seen in Figure 8. As shown in Table 4, all letters for and presented high correlation values and low relative error in quantification and no significant differences between experimental and calculated methods. The voltage maps of letters' distributions shown in Figure 8 (fourth row) are based on voltage intensity distributions, detected by the sensor in the respective positions. In comparison to the voltage maps, the inverse problem reconstructions presented a better spatial resolution, because they solved the letters' geometries, and also allowed a quantitative analysis. Figure 7. Reconstructed distributions of MNP loaded bars' calculated signal. Note that lower voxels were used in the same FOV of 5 × 5 × 1 cm 3 , for a total of 100 voxels of 0.5 × 0.5 × 1 cm 3 . The first row shows the vertical distributions; the second row shows the horizontal distributions; and the third row shows the cross distributions. The columns indicate the distance between the bars' surfaces. The reconstruction of the MNPs' distributions of several cubes arranged to form the letters B, I, O, M, A, and G can be seen in Figure 8. As shown in Table 4, all letters for L exp and L cal presented high correlation values and low relative error in quantification and no significant differences between experimental and calculated methods. The voltage maps of letters' distributions shown in Figure 8 (fourth row) are based on voltage intensity distributions, detected by the sensor in the respective positions. In comparison to the voltage maps, the inverse problem reconstructions presented a better spatial resolution, because they solved the letters' geometries, and also allowed a quantitative analysis. The reconstructed letters' distributions using 361, 100, and 49 scanning points are shown in Figure 9, and the respective quantification parameters are shown in Table 4. The results show that there is not a significant difference between the MNPs' mass quantification using 361, 100, and 49 scanning point reconstructions. The CC values showed a small difference in the quality reconstruction between 361 and 100 points geometries. The 49 points reconstructions showed a considerable reduction in geometry quality in comparison with 361 and 100 points reconstructions. In terms of scanning time, the 361 points scanning took 6 min, the 100 points scanning took 2.5 min, and the 49 points scanning took 1.3 min. Thus, decreasing the number of scanning points may improve the scanning time of ACB, which can improve the method for in vivo applications and the evaluation of physiological dynamic parameters.  The reconstructed letters' distributions using 361, 100, and 49 scanning points are shown in Figure 9, and the respective quantification parameters are shown in Table 4. The results show that there is not a significant difference between the MNPs' mass quantification using 361, 100, and 49 scanning point reconstructions. The CC values showed a small difference in the quality reconstruction between 361 and 100 points geometries. The 49 points reconstructions showed a considerable reduction in geometry quality in comparison with 361 and 100 points reconstructions. In terms of scanning time, the 361 points scanning took 6 min, the 100 points scanning took 2.5 min, and the 49 points scanning took 1.3 min. Thus, decreasing the number of scanning points may improve the scanning time of ACB, which can improve the method for in vivo applications and the evaluation of physiological dynamic parameters.

Discussion and Conclusions
There are currently few biomagnetic techniques capable of reconstructing the spatial distribution of MNPs using non-invasive measurements. Among the techniques reported in the literature, those that have drawn attention and show high potential are MPI, MRX, and MSI. MRX has high sensitivity, can detect micrograms of MNPs, and presents a spatial resolution of a few centimeters [3]. MPI can detect and estimate masses of ferromagnetic material at a microgram scale, with millimeter spatial resolution [35][36][37]. Despite having suitable quality for image reconstruction and quantifying the mass of particles with high sensitivity, both MPI and MRX involve high costs. Thus, the development of alternative biomagnetic techniques capable of reconstructing quantitative MNPs' spatial distributions with low cost, high versatility, and portability, and without the need for electromagnetic shielding, can significantly contribute to scientific research and biomedical applications of MNPs.
Previously, our group showed the effectiveness of the ACB single-channel system in real-time in vivo evaluation of the MNPs' distribution and retention in the liver [19], its circulation in the bloodstream [18], and its different perfusion profile in healthy and injured kidneys [20] in rats. Some of these papers also employed the ACB single-channel system to quantify the MNPs' biodistribution in several organs. Others aimed at the MNPs' cell internalization in ex vivo experiments and also to perform real-time video maps of the

Discussion and Conclusions
There are currently few biomagnetic techniques capable of reconstructing the spatial distribution of MNPs using non-invasive measurements. Among the techniques reported in the literature, those that have drawn attention and show high potential are MPI, MRX, and MSI. MRX has high sensitivity, can detect micrograms of MNPs, and presents a spatial resolution of a few centimeters [3]. MPI can detect and estimate masses of ferromagnetic material at a microgram scale, with millimeter spatial resolution [35][36][37]. Despite having suitable quality for image reconstruction and quantifying the mass of particles with high sensitivity, both MPI and MRX involve high costs. Thus, the development of alternative biomagnetic techniques capable of reconstructing quantitative MNPs' spatial distributions with low cost, high versatility, and portability, and without the need for electromagnetic shielding, can significantly contribute to scientific research and biomedical applications of MNPs.
Previously, our group showed the effectiveness of the ACB single-channel system in real-time in vivo evaluation of the MNPs' distribution and retention in the liver [19], its circulation in the bloodstream [18], and its different perfusion profile in healthy and injured kidneys [20] in rats. Some of these papers also employed the ACB single-channel system to quantify the MNPs' biodistribution in several organs. Others aimed at the MNPs' cell internalization in ex vivo experiments and also to perform real-time video maps of the qualitative MNPs' biodistribution [15,[17][18][19][20]. Furthermore, recently we showed that ACB can assess the effects of corona protein formation in MNPs with three different coatings, which can be employed in biosensing assays [38]. However, the ACB measurements presented to date did not correlate the pixel intensity with the mass of MNPs.
This study presented the ACB quantitative imaging with in vitro and in silico procedures. We also presented a methodology for an ACB scanning approach, which can increase the problem stability and provide quantitative imaging of the MNPs' spatial distribution. By building the ACB sensitivity matrix and using Equation (5), derived from the forward model described in Equation (3), it was possible to estimate the MNPs' mass present in each voxel of the FOV. These procedures may elevate the technique, allowing reconstructions of quantitative images with a better spatial resolution, as shown before for similar MNPs' imaging problems [17].
Using the present scanning approach, we were able to reconstruct the quantitative MNPs spatial distribution in a FOV of 5 × 5 × 1 cm 3 at 5 mm from the sensor surface with 1 cm 3 voxel. In this manner, we obtained a spatial resolution of 1 cm when the SNR is sufficient (MNP masses greater than 1.17 mg), and it was shown that the scanning method has at least 1 cm −1 of resolution and may be improved to reach a higher resolution at a millimetric scale. In this way, the present methodology was able to determine the MNPs' mass and position with good precision, accuracy, and much better spatial resolution (at least 1 cm) than the previous ACB methodologies in the literature, with a relatively high sensitivity [17,21,39]. The measured time for a complete scan was up to approximately 6 min, which can be decreased using other strategies, such as more sensors, continuous scanning, or even by reducing the number of scanning points. In principle, the scanning approach enables an unlimited 2D FOV size. However, increasing the FOV size means increasing the scanning time. It is important to note that as the ACB system is dependent on the MNP's χ(ω), changing the MNPs, or the excitation field frequency, would increase the sensitivity and reconstruct even lower quantities of MNPs. Here, we used a conventional ACB system with a transformer factor of three, and increasing the transformer factor would reflect a higher system's sensitivity [40].
Studies involving MRX have demonstrated the technique's ability to reconstruct tomographic images in a 3 × 6 × 3 cm 3 FOV with 1 cm 3 voxels. MRX can achieve a spatial resolution of millimeters for MNPs' Fe masses higher than a few µg, with an acquisition time of approximately 12 min [26]. In regard to MPI, 3D in vivo scans were performed in a FOV of 20.4 × 12 × 16.8 mm 3 , with submillimeter spatial resolution in horizontal and vertical directions and temporal resolution of a few milliseconds (real-time acquisition). The system showed good sensitivity for tracers in a concentration range of 8 to 45 µmol (Fe)l −1 [36]. It is worth mentioning that for MPI, the voxel size depends on the characteristics and concentration of the MNPs [35]. MSI has shown good results in solving 61 voxels with millimetric volume, and by optimizing the system's features it is possible to solve approximately 108 voxels [39]. Moreover, MSI showed sensitivity up to a distance of 2 cm for milligrams of Fe MNPs [7]. Compared to MRX and MPI, ACB provides similar temporal resolution, but with less sensitivity and spatial resolution and is restricted to 2D reconstructions, at present. However, the ACB system does not demand magnetic shielding and has greater portability. In comparison with MSI, ACB exhibits a similar sensitivity and uses a lower number of voxels for the same sensor to sample distance. Even though MSI and ACB are based on susceptibility measurements, there are substantial instrumental differences between the systems. The ACB is composed of a pair of excitation coils, a gradiometric pair of detection coils, axially arranged, and peripherical electronics optimized for biomedical measurements in vivo [41,42], in vitro [43], and ex vivo [19]. The ACB system's characteristics enable a simpler electronics utilization and also larger FOV imaging. Future works can also focus on an adapted ACB system to quantitatively monitor MNPs in vivo and in real-time [18,19], and perhaps for 3D reconstructions with greater sensitivity and a large number of voxels.
In this work, we presented a new ACB scanning approach that enabled quantitative imaging of in silico phantoms with precisely known MNP masses. The results demonstrated that the scanning method using the ACB mono channel system can obtain quantitative MNPs' images with a spatial resolution of at least 1 cm. The system allows quantifying a few mg of MNPs, which could be further improved by changing the coil's characteristics. In addition, the results obtained through the use of the calculated sensitivity matrix showed that the mathematical-computational model used was adequate, thus allowing new studies based on computational models, given the lower material and temporal cost. These results demonstrate that the ACB system provides quantitative data and can be employed as a tool for MNPs' quantitative imaging of organs and tissues, non-destructively, with a low cost, high versatility, and portability. This data is unprecedented and applicable to other ACB modalities, and will also impact the development of new ACB setups, aiming for in vivo and 3D imaging, which may help the research field relating to MNPs' biomedical applications. Moreover, in the future, it is possible to employ the ACB system and MNPs for biosensing, because it is possible to functionalize the MNPs with specific molecules able to react and link with biomolecules of interest. This can lead to new approaches of ACB to evaluate biochemical bonds of MNPs with plasma proteins and cells, especially using electrochemical-functionalized MNPs for cancer cell targeting, quantitative sensing, and possibly in vivo imaging, including 3D images.