Imaging of Magnetic Nanoparticles With Permeability Tomography

Magnetic nanoparticles are being increasingly used in numerous biomedical applications for diagnosis and therapy. During the course of these applications nanoparticle biodegradation and body clearance may occur. In this context, a portable, non-invasive, non-destructive and contactless imaging device can be relevant to trace the nanoparticle distribution before and after the medical procedure. We present a method for in vivo imaging the nanoparticles based on the magnetic induction technique, and we show how to properly tune it for magnetic permeability tomography, maximizing the permeability selectivity. A tomograph prototype was designed and built to demonstrate the feasibility of the proposed method. It includes data collection, signal processing and image reconstruction. Useful selectivity and resolution are achieved on phantoms and animals, proving that the device can be used to monitor the presence of magnetic nanoparticles without requiring any particular sample preparation. By this way, we show that magnetic permeability tomography may become a powerful technique to assist medical procedures.


Imaging of Magnetic Nanoparticles With Permeability Tomography
A. L. Veiga , M. Fernández-Corazza , M. B. Fernández van Raap , and E. M. Spinelli Abstract-Magnetic nanoparticles are being increasingly used in numerous biomedical applications for diagnosis and therapy.During the course of these applications nanoparticle biodegradation and body clearance may occur.In this context, a portable, non-invasive, nondestructive and contactless imaging device can be relevant to trace the nanoparticle distribution before and after the medical procedure.We present a method for in vivo imaging the nanoparticles based on the magnetic induction technique, and we show how to properly tune it for magnetic permeability tomography, maximizing the permeability selectivity.A tomograph prototype was designed and built to demonstrate the feasibility of the proposed method.It includes data collection, signal processing and image reconstruction.Useful selectivity and resolution are achieved on phantoms and animals, proving that the device can be used to monitor the presence of magnetic nanoparticles without requiring any particular sample preparation.By this way, we show that magnetic permeability tomography may become a powerful technique to assist medical procedures.Index Terms-Magnetic nanoparticle detection, mice tumor, permeability tomography, reconstruction, tomograph design.

I. INTRODUCTION
I N THE last few decades, we have witnessed a growing use of magnetic nanomaterials in scientific research, aimed at biomedical purposes.Nanotechnological approaches try to provide minimal invasive and potentially more effective alternatives than conventional medicine, improving therapeutic outcome and dismissing undesired side effects.In this context, superparamagnetic iron oxide nanoparticles (NPs) have properties that make them valuable tools.Besides displaying low cytotoxicity and biodegradability pathways, they can be externally manipulated with magnetic forces, can release heat under radiofrequency fields, and can be employed as contrast agents.Their large surface/volume ratio allows the functionalization of their surfaces with large variety of molecules, like active targeting ones (antibodies and peptides), drugs, gens, growth factors, and imaging elements like radionuclides and fluorescent molecules [1].Iron oxide NPs have been proposed as carriers of these molecules, cell carriers and heating elements, with promising results in targeted drug delivery [2], localized hyperthermia for cancer treatment [3], [4] and cell therapy for tissue growth and regeneration [5].
All these strategies involve the NP administration systemically or directly at the tissue of interest, followed by its accumulation, distribution and retention at the target tissue during treatment time.Once the NPs are located at the target tissue, repeated treatment may be included in the medical protocol using the already administered magnetic material.But, since the magnetic material may be degraded, for instance in acidic lysosomes [6] and/or be released to body clearance paths [7], the availability of a device that allows magnetic material evaluation during treatments, as well as after treatment completion, is of relevance.
The concern about nanoparticle biodistribution and clearance early led to several propositions for in vivo NPs imaging and evaluation.Among these methodologies is Magnetic Particle Imaging (MPI), first proposed in 2005 [8], that provided the most appealing and useful results.MPI is a tomographic imaging method that guarantees quantification of the spatial distribution of superparamagnetic nanoparticles with high temporal resolution (milliseconds), high spatial resolution (<1 mm), and large sensitivity (μ mol).MPI is based on the relaxation properties of the superparamagnetic nanoparticles.It detects the high harmonics of the magnetization which are strong around a null static field and null at saturation.The technique uses a static field gradient (selection field), an oscillating field (drive field) and signal receive coils.MPI has rapidly advanced from proof of principle to research field and preclinical applications with the first commercially available MPI scanner released in 2014 [9].Even considering the excellent properties and success of MPI, the development of other new imaging techniques, easier to implement, with less expensive instrumentation at the expense of a lower spatial resolution but still useful for research and for a subset of medical problems involving magnetic nanoparticles is desired.
In this work we describe the design, construction and evaluation of a prototype device for imaging magnetic NPs in small animals.The proposed imaging modality is based on Magnetic Induction Tomography (MIT) [10], a technique that uses inductive sensing coils to map the electromagnetic properties of the material under study.
In this technique, an excitation coil is used to generate a primary alternating induction field, producing magnetic polarization and Eddy currents in the material, which induce field perturbations detected by sensing coils.The perturbations can be used for conductivity, permittivity and permeability imaging through reconstruction algorithms.Conductivity imaging has been widely developed for biomedical applications whereas permeability imaging has only been analyzed for industrial applications.
MIT has been analyzed for a broad range of applications as industrial processes, non-destructive examinations and biomedical imaging.The former applications include the inspection of the integrity of steel bars inside of concrete walls and of oil/gas metallic pipes.In biomedical engineering, MIT was implemented by mapping tissue conductivity and applied to hemorrhagic stroke evaluation.A review covering advances of the methodology, instrumentation and applications can be found in [11].One work described the use of a regularized Gauss-Newton scheme for image reconstruction in 3D magnetostatic permeability tomography.The procedure was used to estimate the permeability distribution in phantom-objects built by inserting high and low relative permeability bars (50 and 4) in a cylinder of low relative permeability (3) [14].Another work showed that magnetic permeability of ferromagnetic materials can be reconstructed using frequency difference imaging in the range of 1-100 kHz.Tested objects were carbon steel and electrical steel tubes [15].
To the best of our knowledge, our work is the first design of MIT optimized for the biomedical application of permeability imaging of NPs immersed in living animal tissues.First we show how MIT can be properly tuned to be sensitive to the magnetic NPs, minimizing the contribution of other tissue electromagnetic properties.In the following sections we present the design and construction of the prototype of a tomograph, the reconstruction algorithm, and results on phantoms and in vivo tumor tissues of mice injected with iron oxide NPs, to demonstrate the feasibility and effectiveness of the proposed detection technique.Portability and non-invasiveness were considered mandatory in our design, initially devised to be combined with oncologic hyperthermia treatments [16] to allow the in vivo evaluation of magnetic material persistence at the target tumor.Our solution to evaluate magnetic material in tissues, requires no special sample preparation and provides a fast interpretation of the results that requires no expertise.It is easier to implement than other expensive technologies such as X-ray Computed Tomography (CT) [17], Magnetic Resonance Imaging (MRI) [18], Magnetic Particle Imaging (MPI) [8] and relaxometry measurements RM [19].Among them only RM and MPI are specific for magnetic NPs, being MPI, as already mentioned, the more advanced technique that has lately moved to the preclinical area.
The proposed image generation process is common to other types of tomography, such as electrical impedance tomography (EIT), and consists of two steps: the forward modelling and the reconstruction process.The forward modeling is the simulation of the electromagnetic phenomena based on Maxwell's equations with appropriate boundary constraints.The reconstruction process, also known as inverse problem, is the generation of the tomographic image based on measurements and forward modeling simulations.This is an ill-posed problem and there are many different reconstruction algorithms that address it from different perspectives.Many of them are special cases of the regularized Gauss-Newton approach, including the classic one step linearized Tikhonov regularization with the Laplacian operator [12], [13], [14].

A. Nanoparticles Synthesis and Properties
The magnetic NPs were synthesized according to a previously standardized alkaline coprecipitation method [46], coated with citric acid and resuspended in distilled water to a concentration of 6.3 mgFe/ml.Colloid concentrations were determined by inductively coupled plasma-optical emission spectrometry (ICP-OES, PerkinElmer Optima 2100 DV ICP) at a wavelength of 480 nm [47].The NPs consist of a magnetite core of 13 nm size.Its saturation magnetization is 105 emu/gFe.The citric acid coating endows the NPs with negative charge and ensures electrostatic stability.

B. Agarose Phantom Preparation
This colloid was used to prepare a volume of about 250 μl of agarose phantom 1% (W/V) containing 48 mgFe of NPs in an Eppendorf tube (10 mm diameter).To this end, a gelling mixture was prepared as follows: 1 g of agarose (Invitrogen ultrapure catalogue number 16500-100) was dissolved in 100 ml of phosphate-buffered saline PBS 1X (8 g NaCl, 200 mg KCl, 240 mg KH2PO4 and 1.3 g Na2HPO4 dissolved in one liter of distilled water at physiological pH 7.4) and heated to boiling point until the mixture was translucent.Then, a volume of colloid containing the desired amount of NPs was centrifuged, the NPs were separated using a magnet and resuspended in 75 μl of PBS 1X and mixed in 175 μl of hot gelling mixture and left to cool down to room temperature.

C. Animals
For the in vivo selectivity experiment, mice received two intratumoral injections of 100 μl of the same NPs, suspended in PBS at concentrations adjusted for each animal to incorporate 25 and 50 mgFe of NPs.Two mice received 25 mgFe, labeled as #1 whose tumor size was 11 mm × 10 mm (550 mm 3 ), and #2 whose tumor size was 8 mm × 12 mm (384 mm 3 ).A third mouse, labeled as #3, with a tumor size of 9 mm × 10 mm (405 mm 3 ) received 50 mgFe.The tumour sizes were measured with a caliper and the volumes were estimated as (Dd 2 )/2, where D is the larger diameter and d is the smaller.All experiments with animals followed protocols approved by the Committee of Bioethics and were in accordance with the NIH Guide for the Care and Use of Laboratory Animals.Mice were housed under controlled conditions and food and water were administered ad libitum.Tumor induction is described elsewhere [16].For imaging, the anesthetized mice were placed in supine position in the device and data acquisition was carried out using the same conditions than for the phantom samples.

III. REQUIREMENTS FOR THE DETECTION OF MAGNETIC NANOPARTICLES WITH MIT
Magnetic Induction Tomography is a contactless, noninvasive tomographic technique that aims to map passive electromagnetic properties (PEP) of an object by measuring mutual inductance between pairs of coils placed around its periphery [11].MIT is sensitive to all three PEP properties of the object: electrical conductivity, dielectric permittivity and magnetic susceptibility (σ, and χ m respectively).We show that with a proper control of the design variables, selectivity can be improved for one property of interest, which is magnetic susceptibility in our case.This technique is sometimes named Magnetic Permeability Tomography (MPT) [14]. 1 If proper selectivity is achieved, changes in the susceptibility distribution caused by the presence of NPs should readily be detectable, since magnetic susceptibility of tissues is several orders of magnitude lower than that of therapeutic NP colloids.Magnetic susceptibility of tissues free of NPs is in the order of 10 −6 [26], whereas a therapeutic dose of iron oxide NPs is expected to be in the range from 10 −4 to 10 −3 [27].
A typical magnetic induction array is composed of a number of inductive coils, equally distributed around a circular imaging region, as shown in Fig. 1(a).The procedure consists of exciting one of the coils with an alternate current and sensing the induced voltage in the remaining coils.When an object is placed in the imaging region its properties modify the magnetic field sensed by the non-excited coils.The procedure can be repeated exciting one coil at the time, or rotating the whole array around its center.A tomographic image of a volume can be obtained translating the array along the longitudinal axis perpendicular to the array plane, resulting in a cylindrical imaging region.The calculation of the expected voltage variations in the sensing coils is a complex problem and solutions can be obtained only for simple cases.For example, in [28] an analytical solution was proposed for the case of a cylindrical object of radius R and thickness t positioned coaxially halfway between two coils at distance 2a, with t 2a.We extend this formulation to our case of interest, considering a disk of tissue of radius R, thickness t, and susceptibility close to zero (free of magnetic NPs), including a smaller concentric disk of the same tissue, same thickness but smaller radius R P , which contains non-zero magnetic susceptibility NPs (Fig. 2).Magnetic susceptibility of the tissue with NPs is χ m , being the sole magnetic contribution.The tissue electrical conductivity and relative permittivity are σ and r respectively.The operating frequency is f .Under these conditions, the voltage variations in the sensing coils can be expressed as: where B is the magnetic field, ω is the angular frequency, and K and K P are imposed by the geometry, as: K depends on the dimensions of the whole disk of radius R and K P on that of the disk of radius R P .The thickness t impacts all terms equally, so it has no relative effect on the contributions.The first term of (1) represents the information of interest in our case.It depends on χ m and not on frequency.Second and third terms can be made negligible by selecting a properly low operating frequency, still preserving the first term as a measurable magnitude through advanced instrumentation techniques.This is a particular case of MIT, where the contribution due to induced Eddy currents is intentionally made negligible (null secondary field).The conductivity term is in quadrature and can be separated by phase detection, enabling simultaneous measurement of permeability and conductivity, but as this is of no interest in our case, the quadrature component can be easily eliminated Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.by further reducing the frequency.Frequency reduction has a limit in the current required to produce a measurable magnetic field.Note that reducing the distance between coils is desirable, since K grows faster with a than K P .This distance reduction is limited by the dimensions of the imaging area and is optimized in the final layout.This model is a rough simplification of the more complex geometry of an in vivo situation.However, the analytic formula provides an idea of the order of magnitude of the expected values for the different terms.In particular, we use it to find the conditions under which conductivity and permittivity contributions are negligible and to verify the feasibility of detecting magnetic NPs with our approach.To carry out calculations using equation (1) we assigned to the model depicted in Fig. 2 the dimensions to be used in the final design of our prototype (see details in the next section): R = 40 mm, Rp = 7 mm and 2a = 45 mm.Assigning values of conductivity and permeability for different living tissues, according to [29], [30], in Fig. 3 we show the sensitivity corresponding to the second and third terms of (1) (red and blue lines respectively) as a function of frequency.We found that a susceptibility as low as 0.0001 produces a sensitivity due to the first term (black line) that is 40 dB above the other terms for frequencies below 100 kHz, which constitutes a proper measuring condition.
At the operating frequency, the required sensitivity and selectivity will be verified if the following conditions are met: 1) Signals in the sensing and exciting coils are in phase, evidencing a negligible conductivity contribution.A small copper sample (highly conductive) should produce a visible phase change, while a ferrite sample (highly permeable) should generate an in-phase voltage variation.2) In presence of tissue containing NPs, in-phase voltage increments in the sensing coils must be measurable and must not change its relative magnitude with frequency, evidencing a negligible permittivity contribution.
3) A permeability tomography of a tissue free of NPs must produce an empty image (null voltage variation detected in every sensing coil for any position).In the last section, imaging results obtained with our tomograph prototype on phantoms and living tissues show that these conditions are satisfied.

IV. TOMOGRAPH DESIGN
The analysis shown in the previous section allows us to estimate the required dynamic range of the measurements.Fig. 3 shows that expected values of ΔV /V for coaxial coils are small, in the order of 10 −4 .This situation is aggravated by the fact that the voltage V is much larger in the sensing coils that are located closer to the exciting coil.For example, in a symmetric circular layout as shown in Fig. 1, the voltage induced in coil #2 is more than 10 times the one induced in coil #5 [14].Digital processing of such signals requires an impractically high dynamic range of the analog to digital converters.This is a typical problem in induction tomography, where several solutions have been proposed, like the use of gradiometers [31] or fluxgate magnetometers [32].In order to simplify the layout, we chose to use simple coils combined with a high dynamic reserve lock-in amplifier (> 100 dB).
A second instrumental difficulty arises from the fact that images are reconstructed using the difference between two measurements that can not be registered simultaneously.This is usually called difference imaging or time-difference technique [33], as opposed to frequency-difference methods [15].In our case, induced voltages in the absence of the object are previously measured (V 0 ).Therefore, it is important that long term variations of this voltage due to temperature or other external factors (ΔV 0 ) are not significant compared with variations of V (ΔV ), as induced voltage in the presence of NPs is (V 0 + ΔV 0 ) + ΔV .The measurement is useful only if ΔV 0 ΔV .Therefore, minimization of thermal, mechanical and electrical interaction with the environment becomes relevant in order to increase sensitivity.The long term stability of the system determines how frequently the reference values V 0 must be measured.Hence, a poor stability slows down the imaging procedure.
Our proposal comprises three parts: a magnetically shielded eight-channel coil array, a driver circuit for controlled magnetic field operation, and a lock-in amplifier for raw data acquisition.The diameter of the imaging area (45 mm) is adopted from the dimensions of a portable AC magnetic field applicator we previously designed for oncologic hyperthermia treatment assisted with magnetic NPs in mice [16], [34].This choice provides compatibility among cases to be investigated with both devices.

A. Coil Array and Magnetic Circuit Design
The coils are square shaped and equally spaced in a belt around the imaging area.The number of turns presents a compromise between the required space and the driver power requirements (diameter of the wires).The area of the coils also presents a compromise between the magnitude of the signal (coupling) and the longitudinal resolution (thickness of the slice).In the final design, the coils were designed with 25 turns of isolated  A ferrite core, shown in Fig. 4(a), is used to contain the magnetic flux beyond the imaging region (TDK model B64290L0730X087-N87, outer diameter of 87.0 mm, inner diameter of 54.3 mm, height of 13.5 mm).Its high relative permeability (2200) prevents magnetic interaction with nearby objects, providing a return path for the flux outside the imaging area.A 3D-printed support, shown in Fig. 4(b), was designed for the coils.It is easy to reproduce and, if printed with regular polylactic acid (PLA), has good electromagnetic properties [35].The assembly with the 8 coils mounted is shown in Fig. 4(c).The exciting coil is shielded with a grounded copper sheet to avoid direct capacitive interaction with the sensing coils.Lemo 00 self-latching connectors and coaxial 50 Ω cables are used to avoid interference between channels.

B. Lock-In Amplifier Configuration
As a form of characterization of the array, in Table I we show voltages sensed by the coils when a phantom is placed exactly in the center of the imaging region.Details and discussion on the phantom preparation can be found in the Materials and Methods section.The measurements were taken at 100 kHz, exciting coil #1 (coils numbered as shown in Fig. 1(a) with the maximum 2.5 V sinusoidal output of the lock-in amplifier (Stanford Research Systems Model SR830 100 kHz DSP Lock-In Amplifier).The magnitude, variation and scale in coils #2-8 are presented in the table.The best possible resolution was chosen for each coil in the lock-in amplifier.Note that a small reduction in voltage would have been beneficial, as measuring coils #3 and #7 in the range of 20 mV would produce results with resolution of 1 μV instead of 10 μV.
The second column of the table shows the magnitude of the coupling for each sensing coil.The result is in agreement with previous theoretical analysis in [14].The fourth column shows that relative variations are in the range predicted with the direct calculations we made for a susceptibility of 0.001.A small unbalance can be noted between coils #4 and #6 due to the manufacturing process.A noise voltage of 3.5 mV RMS can be measured in the sensing coils.Given the minimum available lock-in sensitivity of 1 μV, that value represents a signal-to-noise ratio (SNR) of −70 dB.
Placing the phantom exactly in front of any coil produces the maximum coupling between the two adjacent coils.This point of high sensitivity is useful for a fast determination of whether NP selectivity is going to be achieved or not.For example, when driving coil is #1, and a phantom is placed in front of coil #2, a ΔV of 65 μV is measured in coil #3.This value can be used as an upper limit for the measurements of phantoms placed at different position.

C. Driver Circuit and Measuring Procedure
In order to increase its driving capability, the output of the lock-in amplifier was boosted by the circuit shown in Fig. 5.This circuit includes three operational amplifiers connected in parallel to achieve a maximum output current of ±150 mA.A fourth feedback amplifier ensures that the voltage in the excitation coil exactly follows the lock-in reference output.The coil array, Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
the driver circuit and the lock-in amplifier were mounted in a portable disposition to facilitate its relocation to the animal research facility.
The driver circuit is always connected to the same coil, initially in position #1.To obtain a set of measurements, the array is manually rotated, so that the driving coil is placed in positions #1 to #8 in sequence.For each of the positions, the 7 sensing coils are manually switched for reading with the lock-in amplifier.A low resistance (less than 2 Ω) mechanical commutator is used.A total of 56 measurements are taken for each set, which can be represented as a ΔV , 8 by 8, hollow matrix.
Given the noise conditions, the time constant of the lock-in low pass filter must be selected to produce a stable reading.A 300 ms time constant is enough in our case, for which at least 4 s are required for measurements to stabilize.Therefore, the 56 measurements take about 4 minutes to be complete.This time is doubled if a reference slice is taken without the mouse.Some additional time should be considered for translation and rotation of the array.
Several aspects of the design can be modified to increase time resolution.For example, the measurements can be automated replacing the commercial lock-in amplifier by a custom FPGA based design.Also, the rotation of the array can be automated using a simple mechanical configuration, and the manual switching of the coils can be replaced by an electronic analog signal multiplexer.However, the maximum frame rate that can be obtained depends mainly on SNR, which imposes the settling time of the lock-in amplifier.Faster, but lower quality measurements can be obtained by reducing the time constant (increasing noise) and/or by reading a subset of the coils (sacrificing spatial resolution).

V. TOMOGRAPHIC RECONSTRUCTION
The generation of the tomographic images requires the computation of the magnetic field distribution on a virtual model of the device.This computation is governed by Maxwell's equations, but this distribution can be found analytically only for very specific geometrical arrangements.In this section, we describe the methods used to generate the images, starting from a description of the virtual model, followed by the computation of the magnetic field distribution and finishing with the reconstruction algorithm.

A. Virtual 3D Model Construction
We built a tridimensional virtual model of the tomograph accurately following the geometrical dimensions of our real prototype.The tetrahedral mesh consists in three subdomains: the ferrite core (shown in blue in Fig. 6(b), with 102018 tetrahedrons), the air reconstruction domain (translucent gray in Fig. 6(b), 35930 tetrahedrons) and a layer of air (521284 tetrahedrons) filling a rectangular prism of 100x100x40 mm.The mesh was made inhomogeneous to have more resolution in the ferrite core, the reconstruction domain and the air regions near the coils.The volumes of the mesh elements are 0.5±0.3mm 3 for the finer regions and 8±9 mm 3 for the coarser external air layer.The triangular surface boundaries were created using the distmesh package [36] whereas the volumetric tetrahedral mesh was built using the iso2mesh [37] Matlab toolbox.Each of the 8 coils were modelled as a single rectangular coil of zero diameter and 25 turns.

B. Computation of the Forward Problem
The forward problem in MPT is to compute the magnetic field in the space given the geometry, material properties and electromagnetic sources.As previously shown, the magnetostatic approximation holds for the used frequency.This means that the displacement currents and the Eddy current effects can be neglected [11].
We start from the Maxwell equation: In the magnetostatic formulation, the magnetic field H can be decomposed into two fields: the free-space field H s and the magnetization field H m such that H = H s + H m .H s assumes homogeneous permeability and it only depends on the current density sources of the exciting coils J s , whereas H m considers the magnetic properties of the materials and it is irrotational [38]: ( Here we use the reduced scalar potential approach where, due to (5), H m = −∇φ, where φ is a scalar potential.Replacing this identity in (4) we obtain [38]: that must be satisfied everywhere.We apply the far field Neumann condition ∇φ• n = 0 at the outermost boundary of the external air layer, where n is the normal to this boundary vector.Equation ( 6) is equivalent to the volume conduction forward Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
problem in electrostatics, where μ plays the role of the electrical conductivity σ, φ is analogous to the electric potential, and the sources μH s act as electric dipoles in electroencephalography [39] or the external current injection in transcranial electrical stimulation [40] and in EIT [41].The sources in ( 6) are the free-space field H s components that do not depend on the media and can be computed analytically for rectangular coils as we did in this work [45].Then, ( 6) can be solved using conventional finite element method (FEM) code [44].The reduced scalar potential is known to have numerical cancellation errors in the resulting magnetic fields at the high μ domains [38], [43].However, this is not an issue for our particular problem because we are interested in accurate modelling of the magnetic fields at the measurement coils, that are in the low μ air domain.Equation ( 8) with the Neumann conditions can be solved for arbitrary geometries numerically, obtaining the scalar component φ for each mesh node.In this work, we used the classic tetrahedral FEM with linear elements and the Galerkin approach [38], [42].The FEM method converts (8) into a linear system of equations Ka = f , where K is the stiffness matrix accounting for the geometry and the permeability of each layer, a is the unknown potential component φ at each mesh node and f is the independent vector accounting for the electromagnetic sources (in this case, given by H s ).Note that K is sparse and thus, the linear system of equations can be solved efficiently with the preconditioned conjugated gradients (PCG) method implemented in Matlab, using the sparse incomplete LU factorization for the preconditioners.From the computed potential φ, we obtained H m at each tetrahedron by computing the gradient and then B as μ(H s + H m ), recalling that H s was computed analytically.
The measured potential at each measurement coil was computed by following these steps with the aid of the qmeshcut function of iso2mesh: (i) obtain the tetrahedrons that are part of the measurement coil plane and that are inside the coil, (ii) compute the normal to the coil component of B at each of these elements, and (iii) sum them weighted by the area of the section of each tetrahedral element.In other words, we computed the numerical integral of the normal component of B at each coil section, as shown in ( 7) where V i is the measured potential at coil i, N is the number of turns and S i is the section of coil i.Note that (7) represents the forward map from μ to V i which is non-linear with respect to μ due to the H m (μ) = − ∇φ(μ) term.

C. Reconstruction
We applied the one-step classic Tikhonov method to solve the inverse problem [12], which has the general form: In (8), J is the Jacobian matrix, where each column contains the differential simulated measurements obtained by solving the forward problem as described in the previous section for a relative permeability change of 100 at each tetrahedron of a subset of the reconstruction domain.It was arbitrarily selected to have a perturbation in volts of the order of the measurements obtained in the phantom experiments (tens of microvolts).This subset of approximately 8500 tetrahedrons was defined by selecting the elements located at the mid-plane of interest (z = 0) ±2 mm to make computations more efficient.Equation (8) assumes that the problem is linearized around a reference permeability distribution that considers the ferrite core but not the small perturbations produced by the NPs when placed at one tetrahedron at a time.As these perturbations are relatively small, the linearized problem can be considered as a good first approximation, although more iterations might be needed to enhance the precision of the reconstructions [12].Γ is the first-order Laplacian operator in discrete form approximated by finite difference.This means that Γ(i, j) = −1 if i and j are neighbor tetrahedrons (sharing a face), and Γ(i, i) is the number of neighbors of tetrahedron i.This matrix adds spatial smoothness to the solutions and it is the most classical regularization approach [12].α is the hyperparameter, δy are the differences of the measures with and without the NPs, and δμ is the estimated permeability change or tomographic image at each mesh tetrahedron.The hyperparameter was chosen by visual inspection.Note that this is the same algorithm used in a previous MPT work, but using only one iteration [14].Finally, δμ is depicted in a slice at z = 0 of the 3D mesh with no interpolation.

VI. RESULTS
The forward problem was validated by comparing the computed measurements at each of the coils and the real measurements assuming that the phantom is empty.In order to characterise accuracy, we carried out two series of MPT imaging experiments using iron oxide NPs.For accuracy characterization the NPs were uniformly distributed in agarose gels, a material commonly used as tissue equivalent phantoms.Selectivity evaluation was carried out in vivo, imaging NPs injected in mice with tumors.Tomographic parameters and data acquisition were the same for phantoms and for animals.

A. Experiments in Phantoms
In order to optimize the dynamic range, lock-in voltage output was set to a reference peak voltage of 1.7 V for all the experiments.Previous to the measurements, we verified that no signal was detected when samples of living tissue, saline solution or Eppendorf filled with agarose gels were placed in the image region.A syringe loaded with 25 mgFe of NPs was placed in the point of maximum sensitivity and produced a ΔV of 42 μV.
The experiment was conducted placing the Eppendorf tube at three different locations in the internal air area of the device.Fig. 6(a) shows, on the left, the device with a plastic support that has holes to place the Eppendorf tube, which is shown on the right of the same figure.The three locations were (i) the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.geometrical center, (ii) close to a border right in front of coil #2, and (iii) close to a border but in between two coils (coils #1 and #2).The precise locations were measured using a caliper and, considering (x,y) = (0,0) as the central coordinates, the measured center coordinates in mm for the three positions were (0,0), (−8.5, 8.5) and (−4.6, 11.1) [mm], as shown in Fig. 6(c) respectively.
The measurements are shown in the first row of Fig. 7 in the form of matrices where the first row of each matrix corresponds to the measurements obtained when the active coil is #1, the second row are the measurements when the active coil is #2, and so on.Slices at z = 0 of the 3D tomographic images are shown in the last row of Fig. 7 with a white cross marking the peak of the detections.It is observed that the peak of the detections fall inside the phantom tube locations.We conducted the same experiments but using a 2D model, that has the advantage of faster computation times allowing for largest resolution.We found no significant differences between using 2D or 3D models in accuracy, but because some important modelling information might be missed by going from something that is inherently 3D to 2D, we suggest using 3D models whenever is possible.
Note from the measurement matrices in the first row of Fig. 7 that the locations near the borders of the phantom result in a maximum voltage change of 60 μV whereas the location in the center results in a change of 20 μV.In our prototype, the whole array of coils is rotated such that the active coil is physically always the same, instead of switching the active coil as it could be implemented in future versions of the device.We found no differences in the measurements when rotating the coil array with the phantom located at the center, generating a perfect symmetric Toeplitz matrix.This means that the Eppendorf tube was located right at the center of the array, at least considering the precision of our device.The slight misalignment between the reconstruction peak and the geometrical center of the device in Fig. 7 (first row, last column) is due to the coarser mesh used for the 3D model.If interpolation is used, the maximum falls in the node mesh that is closest to the geometrical center.Moreover, in our 2D reconstructions using a much larger resolution, the peak was exactly at the geometrical center.If the active coil is switched (which is not our case), imbalances between different coils could also break the Toeplitz structure of the measurement matrix.Based on the peak position of the reconstructed image, we measured the localization error (LE) as the distance between this position and the central position of the Eppendorf tube.The resulting LE's for positions (i), (ii), and (iii) were 0.7, 1.4 and 0.9 mm respectively which we consider as good enough in terms of accuracy, even when using the most classical reconstruction method.
We conducted a set of simulations to derive resolution metrics and separability of the reconstructed images based on the impulse response or point spread function (PSF), resulting in 7 mm and 12 mm respectively.Sensitivity to iron concentration was also analyzed.Details can be found in Supplementary Material.

B. Experiments in Animals
The in vivo tomographic imaging was performed in order to verify NP selectivity.Before injecting NPs, we verified that the body and the tumor of the mice produced a null matrix.After the injection of NPs, the images of Fig. 8 were obtained.They are intended, previous to ongoing hyperthermia treatment, to verify that the NPs are in place to proceed with the application of the magnetic field.Fig. 8 shows the localization results for the three mice.Our device was able to image the location of the NPs that were injected directly into the tumors of the mice.Most importantly, these results show that the NPs can be imaged when being surrounded by a conductive body, instead of being surrounded by air as in the phantom experiments.Overall, the correct imaging of the NPs in alive animals constitutes the empirical proof that the assumptions previously made are valid and that Magnetic Permeability Tomography is useful for biomedical applications.The image scales and the reconstructed shapes of the domain are indicative, which is typical of this technique, of a better temporal than spatial resolution.The scales are typically considered as unitless, but they could be useful to compare concentrations at different times.It would be even possible to do phantom experiments with different concentrations, sizes and locations of the Eppendorf tubes to calibrate the device and have better estimates of the NP concentrations.Finally, we believe that it would be possible to generate 3D images by moving the object under study or the coil array along the axial direction.Then, each reconstructed image could be stacked to form a 3D image.We plan to add a step motor to perform this movement in the next version of the prototype and study the limits of the technique such as the axial resolution.

VII. DISCUSSIONS
The proposed MPT is a technique that exclusively maps magnetic permeability without including any contrast agent.Since the permeability of therapeutic nanoparticles is several orders of magnitude larger than that of tissues, they themselves work as contrast.The images the method produces just show the nanoparticles with minimal interaction with other properties of the tissues.
The basic principle of MPT is different from other methodologies also proposed for tracing NPs, which are quite appealing techniques, but require a more complex, expensive, and bulky instrumentation.For example, Magnetic Relaxometry uses Superconducting Quantum Interference Device (SQUID) sensors Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
to measure NP magnetic relaxation.MPI uses the nonlinear magnetization response of NPs to an alternating magnetic field while a field gradient is scanned.CT has been used to visualize deposits of nanoparticles in prostate tissue during magnetic hyperthermia clinical applications while MRI could not be used due to signal loss in the regions containing the NPs [21].More recently, sweep imaging with Fourier transformation (SWIFT) MRI technique showed a linear relation between longitudinal relaxation and iron concentration up to 3 mg Fe/ml [18] and MR susceptibility maps has been reconstructed from MR phase data and applied to detect magnetic NPs [20].These are quite appealing techniques, but require a more complex and expensive instrumentation.Magneto acoustic tomography (MAT) has also been proposed as an imaging alternative to MRI or CT to image nanoparticle distribution in tissues.The application of magnetic pulses, of microsecond duration, generates pulsed magnetomotive forces that act on the NPs creating acoustic vibrations that propagate in all directions of the medium.The ultrasound signals are registered in a circular configuration, and tomographic images are reconstructed [22], [23], [24].EIT uses a similar experimental set up than permeability tomography, but is meant to measure electrical conductivity differences across sections of the organs or tissues.Moreover, it requires the modeling of the tissues because the conductivity of the nanoparticles is in the same order of the surrounding tissues.In contrast, magnetic permeability of the nanoparticles is several orders stronger than the magnetic permeability of biological tissues, which constitutes a major advantage of permeability tomography against EIT.
The main drawback of MPT is related to the solution of an inverse problem to produce the image after data collection.Image spatial resolution can be poor, and the images look blurry.But as can be seen in Fig. 8, due to the high permeability difference between NPs and tissues, the eight coil configuration is enough to detect where the NPs are accumulated.The accuracy estimated from the phantom experiments and the spatial resolution estimated with the simulations are enough to verify the feasibility of the proposed technique and its purpose: to confirm the presence and location of the implanted nanoparticles inside the tumor to be treated with a simple, portable and low-cost method.This resolution can be improved using a set up with more coils.For more precise temporal surveillances, calibration with known density nanoparticle samples is required.Also, blurriness and precision can be improved to some extent with the use of more advanced reconstruction algorithms such as total variation [12], [13] or spatial filters [41], [48].However, comparison of different reconstruction algorithms is out of the scope of this work where we want to highlight that even the simplest method is good enough as proof of feasibility for this application.Nonetheless, MPT displays significant advantages for medical use.The data acquisition process is non-invasive and contactless.The instrumentation is relatively easy to implement and not expensive.No ionizing radiation is used with MPT (compared with CT) and the magnetic material is sensed without any labelling as it happens with single photon emission tomography SPECT and positron emission tomography PET, which can be used to analyze magnetic material distribution provided that the NPs are labeled with an appropriate radionuclide [25].There are no risks related to their use as happens with MRI strong magnets.The tomograph is light and compact enough such that it can be easily transported.
One limitation of this study might be the use of a one-step reconstruction algorithm and the assumption that the linearization of the forward problem is a good approximation for the small perturbation caused by the nanoparticles.However, the phantom experiments showed that the results were accurate enough for this first proof-of-concept of the technique, even when using a very simple reconstruction algorithm.We plan a future study analyzing the performance of using more steps in the reconstruction algorithm as well as other reconstruction methods such as total variation [12], [13] or spatial filters [41], [48].In future works, phantoms with different NP concentrations can be used to analyze in detail the performance of different reconstruction algorithms that might improve accuracy, spatial resolution and shape of the reconstruction.

VIII. CONCLUSION
We presented a proof of concept of an MIT device for imaging NPs in a living tissue.To this end, we first showed theoretically how device parameters must be tuned to increase permeability selectivity and make it less sensitive to conductivity and permittivity.We presented our tomograph prototype and provided exhaustive technical construction details and a detailed discussion on imaging reconstructions procedures.The proposed imaging method was successfully tested in phantom and in vivo experiments, showing that it may constitute a powerful technique to monitor NPs incorporation in the target tissue and its persistence in applications like magnetic hyperthermia treatments.To our knowledge, this is the first MIT prototype specifically designed for this biomedical application.

Fig. 1 .
Fig. 1.(a) Eight coil array with an object of interest in the imaging region.(b) Static simulation used to estimate flux density redistribution (Vizimag visualization software version 3.185) when coil #1 is excited and a high permeability shield is placed in the periphery.(c) Tomographic reconstruction from sensed voltages when AC voltage is applied to each coil in sequence.

Fig. 2 .
Fig. 2. Simplified model used for a direct estimation of the induced voltage.A tissue disk of radius R, including a concentric region of radius R P , that represents a tissue loaded with NPs,is placed in the middle of two coils separated a distance 2a.Thickness of the disk t is much smaller than a.

Fig. 3 .
Fig. 3. Contributions to relative voltage variations of conductivity (red lines) and permittivity (blue lines) of different tissues such as blood, liver, heart, fat, kidney, lung, and intestine.At 100 kHz, a susceptibility value of 0.0001 (black line) still provides a useful 40 dB selectivity margin.

Fig. 5 .
Fig. 5. Driver circuit for the exciting coil.Resistors of 3.3 Ω are included to equalize the outputs of the amplifiers.The 5.6 kΩ/3.3 pF network is included to achieve stability.Powered by a ±3 V power-supply, this circuit operates properly at 100 kHz.It delivers a current up to ±150 mA, which is enough to drive a 23 μHy coil with a 2 V peak sinusoidal voltage.

Fig. 6 .
Fig. 6.(a) Coil array with connectors and slots for the Eppendorf tubes and sample holder is shown on the right.(b) Virtual 3D model of the device where the ferrite core is shown in blue and the domain of reconstruction in translucent gray.The rectangular coils are depicted with a red line and numbered in (c) where also the positions used for the phantom experiments are shown as red circles.

Fig. 7 .
Fig. 7. Results of phantom imaging experiments for positions (i), (ii) and (iii) of Fig. 6(c).The upper panel shows the matrices of the voltage measurements [μV] where the matrix row numbers indicate the active coil and the column numbers the measurement coils.The lower panel shows the reconstructed images, where the white crosses indicate the peak location and the red circles the ground truth locations, same as in Fig. 6(c).

Fig. 8 .
Fig. 8. Results of in vivo imaging experiments for three different mice.The upper panel shows pictures of the animals where approximate locations of the tumors are highlighted with white lines.The middle panel shows the matrices of the voltage measurements [μV] where the matrix row numbers indicate the active coil and the column numbers the measurement coils.The lower panel shows the reconstructed images.The first two columns of graphs correspond to experiments with 25 mg of NPs and the last column with 50 mg of NPs.The scales of the first two columns differ from the third one.

TABLE I MEASUREMENT
RESULTS OF A PHANTOM IN THE CENTER OF THE IMAGING REGION, WHEN EXCITING COIL #1 AT 100 KHZ.FOR EACH COIL, EITHER SCALE A (20 MV RANGE, 1 μV RESOLUTION) OR B (200 MV RANGE, 10 μV RESOLUTION) WAS SELECTED IN THE LOCK-IN AMPLIFIER TO BE AS CLOSE AS POSSIBLE TO FULL-SCALE OPERATION 0.2 mm copper wire in two layers, 12.5 mm × 15.8 mm area and 2.6 mm height.Inductance of the coils is 23 μHy and its resistance is 0.9 Ω.In order to characterize manufacturing repeatability, each coil was measured at 100 kHz, resulting in a maximum inductance value of 23.85 μHy and a minimum of 22.59 μHy.