Human-sized quantitative imaging of magnetic nanoparticles with nonlinear magnetorelaxometry

Objective. Magnetorelaxomety imaging (MRXI) is a noninvasive imaging technique for quantitative detection of magnetic nanoparticles (MNPs). The qualitative and quantitative knowledge of the MNP distribution inside the body is a prerequisite for a number of arising biomedical applications, such as magnetic drug targeting and magnetic hyperthermia therapy. It was shown throughout numerous studies that MRXI is able to successfully localize and quantify MNP ensembles in volumes up to the size of a human head. However, deeper regions that lie far from the excitation coils and the magnetic sensors are harder to reconstruct due to the weaker signals from the MNPs in these areas. On the one hand, stronger magnetic fields need to be applied to produce measurable signals from such MNP distributions to further upscale MRXI, on the other hand, this invalidates the assumption of a linear relation between applied magnetic field and particle magnetization in the current MRXI forward model which is required for the imaging procedure. Approach. We tackle this problem by introducing a nonlinear MRXI forward model that is also valid for strong magnetic excitation fields. Main results. We demonstrate in our experimental feasibility study that scaling up the imaging region to the size of a human torso using nonlinear MRXI is possible. Despite the extreme simplicity of the imaging setup applied in this study, an immobilized MNP sample with 6.3 cm3 and 12 mg Fe could be localized and quantified with an acceptable quality. Significance. A well-engineered MRXI setup could provide much better imaging qualities in shorter data acquisition times, making nonlinear MRXI a viable option for the supervision of MNP related therapies in all regions of the human body, specifically magnetic hyperthermia.


Introduction
Magnetorelaxometry imaging is a noninvasive and nonionizing experimental imaging procedure for localizing and quantifying magnetic nanoparticles inside a volumetric region of interest (ROI) (Flynn and Bryant 2005, Liebl et al 2015. Alongside magnetic particle imaging (MPI) (Gleich and Weizenecker 2005) and magnetic resonance imaging (MRI) (Paysen et al 2019), Magnetorelaxomety imaging (MRXI) presents another interesting approach to monitor magnetic nanoparticle (MNP) related therapies, such as magnetic drug targeting (Alexiou et al 2011, Lyer et al 2015 and magnetic hyperthermia therapy (Thiesen and Jordan 2008, Bañobre-López et al 2013, Salunkhe et al 2014. The in vivo location of the MNPs obviously plays a major role in the efficacy of the respective therapeutic approaches. Additionally, the MNP quantities are also of importance since they directly influence the drug concentration in magnetic drug targeting and the heat dissipation during magnetic hyperthermia therapy. MRXI consists of two essential steps. First, the MNPs inside the ROI are aligned along a static magnetic field generated by externally arranged electromagnetic coils to induce a net magnetization. In the second step, the coils are rapidly deenergized, switching off the magnetization field so that the magnetic moments of the MNP ensemble can fall back into its equilibrium state by thermal agitation which results in a decay of the net magnetization over time. This entire physical process is mathematically described by the MRXI forward model (Baumgarten et al 2008, Liebl et al 2014 and by measuring the magnetic (stray) field of the decaying net magnetization (i.e. the relaxation), it becomes possible to reconstruct the MNP ensemble quantitatively by solving the corresponding ill-posed inverse problem. It has been shown that MRXI is capable of localizing and quantifying MNPs in dm 3 -sized volumes (Liebl et al 2015, Arsalani et al 2022. Nonetheless, deeper lying regions of a large ROI are difficult to reconstruct due to the very weak magnetic excitation fields that are typically applied in MRXI (Arsalani et al 2022). From a technological point of view, further upscaling of the ROI would merely require stronger magnetic excitation fields and suitable sensors to record the smaller relaxation signals of more distant sources. These fields could still easily be weak enough to not trigger arising health risks such as tissue heating or peripheral nerve stimulation in the patient (Panagiotopoulos et al 2015, Billings et al 2021. However, due to the resulting large magnetic fields that drives MNPs into their nonlinear magnetization regime (Sarangi et al 2009), one initially needs to adapt the state-of-the-art MRXI forward model that is designed for linear magnetizations.
The goal of this feasibility study is to demonstrate that MRXI can be scaled up to the size of a human torso with only few hardware requirements. To our knowledge, this is the first study that attempts to image MNPs in such a large volume using MRXI. We aim to achieve this by using a large excitation coil that is able to appropriately magnetize a torso-sized volume. For this purpose, we introduce the nonlinear MRXI forward model that is able to operate in a nonlinear MNP magnetization regime. We use an optically pumped magnetic gradiometer (OMG) for our MRXI experiment. This type of sensor is much more flexibly positionable than the commonly employed SQUID-systems (Crevecoeur et al 2012, Liebl et al 2014, that are spatially constrained by their required liquid helium cooling system, and is thus conveniently usable out of the box in our experiment. By providing this proof of concept we aim to verify that MRXI is a reasonable option to monitor MNP related therapies in large human body parts.

Components of the experimental human torso MRXI setup
The entire measurement setup consists of a single OMG, a single excitation coil, an MNP sample, as well as a Helmholtz coil pair. Furthermore, it is located inside a magnetically shielded room with a shielding factor of approximately 100 at 0.01 Hz and exceeding 10 000 above 10 Hz with inner spatial extensions of about 1.5 m × 2 m × 2 m. A photograph of the measurement setup from inside the magnetically shielded room is shown in figure 1.
The OMG we use was developed by Twinleaf, USA. We employ a gradiometer rather than a single-channel optically pumped magnetometer since the residual background field can be filtered out much more efficiently in this arrangement. The two optically pumped magnetometer cells inside the OMG are 2.3 cm apart and the front cell lies only 5.6 mm behind the front surface of the housing, meaning the sensor can be positioned very closely to a patient. The measurement sensitivity in gradiometric mode lies slightly below 0.2 pT cm Hz ( )and the sampling frequency is set to 1 kHz. The clock signal for the sampling frequency is generated by an AFG-2225 Dual-channel arbitrary function generator, GW Instek, Netherlands, which is also used to trigger the magnetization and relaxation phases during the MRXI experiment with its second channel (see section 2.4). Additionally, a Helmholtz coil pair with a radius of 0.7 m is operated inside the magnetically shielded room, generating a static and homogeneous magnetic field of approximately 15 μT/μ 0 in its center point. This is necessary because the OMG requires a stable background field along its sensitive direction in order to function properly. Such moderate magnetic fields are assumed to be unproblematic for the generation of the MNP relaxation, just causing an offset in the measured magnetic fields during MRXI (Jaufenthaler et al 2020).
The excitation coil consists of 160 windings with an inner diameter of 5.6 cm and an outer diameter of 25.8 cm. At the maximum drive current of 60 A that is applied during our experiments the magnetic field strength reaches approximately 70 mT/μ 0 in the coil center point. The rise time of the coil current is approximately 28 ms, after which 90% of the peak amplitude are reached. The fall time is much shorter due to transient voltage suppression diodes that are connected in parallel to the coil. It lies around 3 ms, after which the field has dropped to 10% of the peak amplitude.
It was shown in in vitro and in vivo experiments that the magnetization decay of the MNPs located inside organs or tumor tissue is governed by the Néel relaxation process (Richter et al 2010, Dutz et al 2011, Farzin et al 2020. Therefore, the sample we use for our experiments are 12 mg Fe of perimag ® (plain, micromod, Germany) nanoparticles immobilized in an ellipsoidal block of gypsum sized 3 cm × 2 cm × 2 cm in the three spatial dimensions. This dosage equates 0.08 mL of MNPs per 1 mL tumor volume which relates to a very low amount for effective magnetic hyperthermia therapy (Johannsen et al 2005(Johannsen et al , 2010. Since dosages administered in clinical practice are typically much higher and the magnetorelaxometry signal intensity is proportional to the MNP iron mass, this sample presents a worst case scenario for MRXI.

Nonlinear MRXI model
The linear forward model has been commonly applied for MRXI throughout recent years and is described in detail elsewhere (Crevecoeur et al 2012, Liebl et al 2014. A brief overview is given here to be able to follow the conducted adaptations that are necessary to form the nonlinear MRXI forward model.
A typical MRXI setup consists of multiple excitation coils, N s magnetic sensors and a region of interest that is tessellated into N v volume elements or voxels. An MNP ensemble within the ROI is considered via a certain iron mass c v per voxel v. After magnetizing the MNP sample for an appropriate time window t mag using the excitation coils, the coils are deactivated abruptly and the magnetization decay is recorded by the N s sensors. The relaxation amplitudes B s between two arbitrarily selected time points t 1 and t 2 after the deactivation of the coils are then extracted from the measurements and can be modeled for each sensor s with where μ 0 denotes the magnetic permeability of free space, n s the unit vector describing the sensitive axis of sensor s, and r s,v the vector pointing form voxel v to sensor s. The factor κ 1,2 = κ(t 1 ) − κ(t 2 ) is the amplitude decay between the time points t 1 and t 2 of the relaxation function κ(t) that phenomenologically describes the MNP relaxation and is normalized to a range between 0 and 1. κ(t) can be modeled logarithmically (Chantrell et al 1983), as stretched exponential (Eberbeck et al 2006), or as a sum of exponential functions  as done in this study. Lastly, t mag ( ) c denotes the magnetization time dependent magnetic susceptibility of the MNPs and h v is the magnetic field vector in the center of voxel v that is produced by the excitation coils.
The factor t mag 1,2 ( ) c k can be acquired through characterization measurements of the administered MNPs and is approximated by a constant value in the linear MRXI forward model . Due to the stronger magnetic excitation fields that are necessary to sufficiently magnetize MNP samples inside an ROI as large as a human torso, this approximation is no longer valid since some areas are inevitably driven into the MNPs' nonlinear magnetization regime. This affects the relaxation times (Sarangi et al 2009, Deissler et al 2014 as well as the magnetization (Bedanta and Kleemann 2008) of the MNPs and consequently influences both the magnetic susceptibility χ and the phenomenological relaxation behavior κ. Therefore, this magnetic field dependency must be accounted for in the factor and thus gives rise to the nonlinear MRXI forward model Ît he system matrix. In order to simulate the gradiometric measurements as provided by the OMG, L is modeled as the difference of the system matrices of the front and the rear sensor cells with 2.3 cm baseline in between (see section 2.1) such that L = L front − L rear .
2.3. Recording the magnetic field dependent factor As indicated in section 2.2, the magnetic behavior of the MNPs must be characterized to enable quantitative imaging using MRXI. Specifically, this means that the magnetic field dependent factor f h v 2 ( )   has to be recorded over the entire range of magnetic field amplitudes h v 2   applied in the ROI of our imaging setup (10 μT/μ 0 up to 30 mT/μ 0 ).
The characterization setup is shown in figure 2. The sample is placed centrally in front of the OMG at a distance of 7 cm from the face of the OMG to the ellipsoid's center. The excitation coil is placed centrally on the opposite side of the MNP sample with distances of 3 cm, 8 cm, 13 cm, 18 cm, and 28 cm throughout 5 individual measurements to appropriately cover the entire range of the magnetic field amplitudes of interest and simultaneously to assess the validity of the simulation model with respect to the calculated field amplitudes. The coils are driven with coil currents between 4 A and 60 A within 7 linearly spaced steps. The magnetization time is t mag = 200 ms and the subsequent relaxation time is t relax = 300 ms for every excitation, resulting in a 500 ms cycle time for each relaxation amplitude measurement. Additionally, empty measurements without MNP inside the chamber are recorded for every coil position and subtracted from the MNP relaxation measurements to eliminate the influences of residual magnetic fields from the measurement environment. Every relaxation amplitude measurement is averaged over 8 repetitions for each coil current and each coil position.
Analogous to , the magnetic field dependent factors f h i 2 ( )   per relaxation measurement i are then calculated as ratios between measured relaxation amplitudes B i,meas and the corresponding simulated relaxation amplitudes B i,sim (see also (3)) such that

Imaging setup and procedure
The imaging setup we apply is shown in figure 3. It was merely designed to show the feasibility of imaging a torso-sized volume (slightly over 20 cm × 30 cm × 40 cm) using MRXI with only a few components. It is also located inside the magnetically shielded room and consists of a single OMG and a single excitation coil, coaxially placed at a distance of 25 cm between the front faces of coil and OMG. The MNP sample holder is connected to the torso baseplate such that the MNP sample can be placed approximately at the position where a pancreatic tumor in a patient might be located. We chose a pancreatic tumor since its locally advanced and borderline resectable variants present interesting targets for magnetic hyperthermia therapy (Attaluri et al 2020). A lifesized torso volume was 3D-printed and fixated on the torso baseplate, representing the considered ROI for the imaging experiment.
Commonly, an MRXI setup consists of a multitude of sensors and excitation coils. We mimic such a typical MRXI setup by moving the torso to various spatial positions, each time measuring the relaxation amplitudes of the MNPs using the same excitation procedure as described in section 2.3. The torso baseplate is moved throughout 4 positions in 9 cm steps in y-direction on the grid beneath, rotated 180°around its central z-axis, and moved throughout the same 4 y-positions again. This procedure is repeated on 4 different z-positions where the baseplate is incrementally raised 5.4 cm each time. This results in 32 different position pairs of the coil and the OMG around the ROI. Figure 4 depicts the MRXI setup modeled in a MATLAB simulation environment. The ROI is a 20 cm × 30 cm × 40 cm sized cylinder that is discretized into 8 × 12 × 16 voxels (gray cubes) with an edge length of 2.5 cm each. The 32 coil positions that result from moving the torso through space are shown in different colors ranging from blue to yellow. The corresponding OMG positions are depicted as red arrows on the opposite side of the ROI, respectively. The coil is driven by 7 different coil currents at every position, resulting in a total of  b meas 224 Î measured relaxation amplitudes. Since every voxel holds a specific MNP iron mass c v , the vector we seek to reconstruct is of the size  c 1536 Î , making our system of linear equations (4) heavily underdetermined and accordingly hard to solve.
We use an elastic net regularization (Zou and Hastie 2005) with additional non-negativity constraint to solve the inverse problem of (4). The elastic net regularization, which combines L 1 -and L 2 -penalization, is a reasonable choice for MRXI (especially for an ROI as large as a torso) since the MNPs are typically concentrated within relatively small subvolumes which results in a sparse MNP mass vector c. Also, the variables in MRXI are often highly correlated where the elastic net performs better (Zou andHastie 2005, Bühlmann et al 2013) than standard regularization methods, such as Tikhonov regularization (Golub et al 1999) or LASSO (Tibshirani 1996).

Recording the nonlinear magnetic field dependency
The characterization measurements described in section 2.3 yield a total of 35 averaged MNP relaxation measurements. The relaxation amplitudes are extracted as the signal difference between the time points

Imaging results
The total measurement time for the procedure described in section 2.4 is approximately 15 min including the empty measurement. Naturally, the relaxation amplitudes of the 224 averaged pancreatic tumor relaxation measurements are also extracted using the time points t 1 = 40 ms and t 2 = 190 ms to ensure the same excitation field dependency as shown in figure 5(b). The measured relaxation amplitudes range from a few pT up to a few hundred pT. This yields a signal-to-noise ratio of only SNR P b 10 log 13 dB, 6 meas 2 2 noise · ( )   = » since the sensor noise is of similar magnitude (on average about 30 pT peak-to-peak). P noise denotes the average power of the sensor noise. Nonetheless, it is possible to localize and reconstruct the MNP distribution inside the tumor volume in a decent quality as depicted in figure 6. The left images show the ground truth MNP distribution as used in the imaging experiment in different orientations, respectively. High MNP masses per voxel are shown in dark red and low MNP masses in light red. Additionally, the voxels become increasingly transparent with decreasing MNP mass. The ellipsoidal shape of the sample is lost to the coarse discretization of the ROI. The right images show the reconstructed MNP distribution yielded by the elastic net regularization with the nonlinear MRXI forward model and the measured relaxation amplitudes as inputs.
We assess the qualitative similarity between the reconstructed and the ground truth MNP distribution via the Pearson correlation coefficient CC (Starovoytov et al 2020), where CC = 0 indicates no similarity at all and CC = 1 means a perfect match. Using our very simple MRXI setup, a qualitative overlap of CC = 0.81 can be reached for the pancreatic tumor MNP localization. In total, an MNP mass of 13.1 mg is reconstructed, which equates a relative deviation of 5.2% from the ground truth phantom mass.
In contrast, using the typical linear MRXI forward model with a constant magnetic field dependent factor of f h 4.8 10 m kg v 2 4 3 1 ( ) ·   = -yields a dramatically decreased model fidelity. As a consequence, it is not possible to localize the MNP ensemble (CC ≈ 0) and the quantification deviation exceeds 70%.

Discussion
Throughout this study, we show that the reconstruction of an MNP ensemble within a human torso-sized volume is feasible for nonlinear MRXI, already with a rather crude setup design as the one applied here. The general localization of the sample works well and also the quantification of the MNPs yields a good result with only 5.2% deviation from the correct value. The incorporation of the nonlinearity into the MRXI forward model is crucial for large magnetic excitation fields as proven by the ineligible reconstruction yielded by the linear MRXI forward model.
However, even with nonlinear MRXI there remain inaccuracies in reconstructing the exact shape of the sample, which were expected given the simplicity of the imaging setup and the severely underdetermined forward model. Nonetheless, despite the crudeness of the setup, a reasonable correlation coefficient of CC = 0.81 was reached. The remaining inaccuracies can be strongly mitigated by using multiple sensors as in most other MRXI studies (e.g. (Flynn and Bryant 2005, Adolphi et al 2012, Liebl et al 2014 to name a few), and by designing optimized excitation coil grids  as well as optimized excitation strategies (Schier et al 2019). Also, using Superconducting Qantum Interference Devices (SQUIDs) instead of OMGs would contribute considerably to raising the SNR above merely 13 dB due to their lower noise levels. The OMG has a sensitivity of approximately 0.2pT/(cm Hz ) while modern SQUID systems operate at noise levels around 1 fT/ Hz (Schnabel et al 2004, Storm et al 2016. On the other hand, OMGs are much more flexibly placeable which would necessitate finding a good trade-off between sensitivity and positioning for better MRXI setups.
A voxel edge length of 2.5 cm is the lower limit with the current setup that reliably yields a decent reconstruction quality of the MNP distribution. Higher ROI resolutions would also require a more sophisticated MRXI setup. Specifically, because only one coil is applied in the current setup, it was designed to provide a good magnetic field coverage of the entire ROI and is accordingly relatively large. A grid of smaller, more densely wound coils would produce more inhomogeneous magnetic fields that decrease the linear near-dependency of voxels close to each other and thus make them easier separable in the solution of the inverse problem. Using this approach, we believe that the resolution of the ROI can be increased considerably.
The measurement time for the current imaging procedure is about 15 min. This rather long duration can be decreased substantially by reducing the number of different coil currents, measurement averages, as well as the t mag and t relax times. The 7 different coil currents are mainly applied to gather additional data points for the reconstruction from the nonlinear behavior of the system. However, the relaxation responses differ only slightly from each other and thus contribute marginally to the final reconstruction quality. It would be much more beneficial for both the measurement time and the reconstruction quality to simply use more (at best optimized (Schier et al 2022)) coil positions, each of which individually driven by only one or two coil currents, or by using optimized excitation strategies (Schier et al 2019). The necessity of averaging relaxation amplitudes can be largely attributed to the repetition inconsistency of the coil current supply hardware. The coils are driven by a very simple laboratory power supply which is not built for inductive loads. Through repeated relaxation measurements it was empirically determined that averages of 8 measurement trials yield an adequate repetition accuracy. However, a coil current supply that is specifically designed for driving coils with rectangular pulses might even eliminate the need for trial repetitions, thus vastly decreasing the overall measurement time. Finally, the durations t mag and t relax are not optimized but arbitrarily chosen. It is safe to assume that these time windows can be shortened to some extent as well.
A minor uncertainty regarding the MNP quantification in a realistic imaging scenario still remains after our experiment. Since we used the same sample for characterization of the magnetic field dependent factor f h i 2 ( )   as for the imaging experiment, only a small deviation between the reconstructed and the nominal MNP mass was to be expected. In reality, one typically does not have access to an MNP sample with binding state properties identical to the ones prevailing in the (pancreatic) tumor of the patient. However, since the majority of MNPs in tumor tissue are immobilized (Richter et al 2010, Dutz et al 2011, Farzin et al 2020 and thus governed by Néel relaxation which is mostly unaffected by environmental parameters, we do not expect a considerable deterioration of the viable MNP quantification as long as the magnetic characterization is performed using identical, immobilized MNPs at the appropriate (body) temperature. Even if such deviations in the relaxation would be present, they could be identified in the shape of the MRXI curves and be corrected by using a more appropriate sample for determining f h i 2 ( )   .

Conclusion
We successfully demonstrated imaging of a torso-sized volume by nonlinear MRXI is superior to linear MRXI at higher excitation fields. Although the imaging setup applied here is far from optimal, it enables localization and quantification of an immobilized MNP sample with 6.3 cm 3 and 12 mg Fe with an acceptable quality. A wellengineered MRXI setup could provide much better imaging qualities in shorter data acquisition times, making nonlinear MRXI a viable option for the supervision of MNP related therapies in all regions of the human body. Specifically, magnetic hyperthermia might benefit tremendously from such an imaging procedure since its spatially slowly changing MNP distributions are especially well suited for monitoring by MRXI.

Data availability statement
The data cannot be made publicly available upon publication because no suitable repository exists for hosting data in this field of study. The data that support the findings of this study are available upon reasonable request from the authors.

Funding
Financial support by the Austrian Science Fund (FWF, grant I 4357-B) and the German Research Foundation (grant WI 4230/4-1) is gratefully acknowledged.