Evaluating selection criteria for optimized excitation coils in magnetorelaxometry imaging

Objective. Magnetorelaxometry imaging (MRXI) is an experimental imaging technique applicable for noninvasive, qualitative and quantitative imaging of magnetic nanoparticles (MNPs). Accurate reconstructions of nanoparticle distributions are crucial for several novel treatment methods employing MNPs such as magnetic drug targeting or magnetic hyperthermia therapy. Hence, it is desirable to design MRXI setups such that the reconstruction accuracy is maximized for a given set of design parameters. Several attempts exist in literature that focus on the improvement of MRXI and other related linear inverse problems with respect to various figures of merit. However, to date it remains unclear, which approach leads to the largest benefit for the reconstruction accuracy. Thus, the aim of this study is to compare the different figures of merit, thereby determining the most reliable and effective optimization approach for magnetorelaxometry setups. Approach. In the present simulation study, we translate these figures of merit to various cost functions, allowing us to optimize the electromagnetic coil positions and radii of two distinct MRXI setups with an adapted tabu search algorithm. Multiple artificial MNP phantoms are reconstructed employing the optimized setups and the resulting imaging qualities are subsequently compared. Main results. The extensive amount of generated synthetic data unprecedented in previous MRXI studies identifies the condition number as the most reliable indicator for good imaging results. This is the case for both the qualitative as well as the quantitative reconstruction accuracies. Significance. The results of this study show that optimized coil configurations increase the reconstruction quality compared to the state-of-the-art. The insights obtained here can also be extended to other design parameters of MRXI setups, thus enabling more reliable reconstructions of MNP ensembles which will ultimately render the aforementioned treatment methods safer and more efficient.


Introduction
Magnetic nanoparticles (MNPs) are an emerging class of nanomaterials that are widely applicable for novel biomedical therapies and treatment methods due to their unique properties. In particular, these approaches include magnetic hyperthermia (Thiesen andJordan 2008, Bañobre-López et al 2013), magnetic drug targeting and controlled drug release (Alexiou et al 2011, Lyer et al 2015, tissue engineering (Ito and Kamihira 2011), as well as lab-on-a-chip technologies (Medina-Sánchez et al 2012). Especially the field of oncology might benefit tremendously from the novel options in tumor therapy. An excellent overview about the recent advances of MNP employment in biomedicine can be gained from Cardoso et al (2018). However, a prerequisite for the success of in-vivo MNP applications such as hyperthermia therapy and drug targeting is to have precise knowledge about the location of the particles for these approaches. Additionally, it is required to be able to quantify the MNP concentration in these areas since it is directly related to the heat dissipation during hyperthermia therapy and the amount of medication released during drug targeting. Hence, the ability to monitor these parameters ensures the safety and efficiency of the aforementioned treatments.
Magnetorelaxometry imaging (MRXI) is an experimental imaging technique that is able to accurately reconstruct MNP distributions qualitatively as well as quantitatively in a noninvasive manner , Liebl et al 2015. MRXI utilizes the concept of magnetorelaxometry (MRX) (Kötitz et al 1995) in which the time-delayed response of the superparamagnetic nanoparticles (i.e. the relaxation signal) to several previously applied, spatially different magnetic field configurations is measured. Under certain simplifications, the entire process of MRXI can be modeled as a linear system of equations, referred to as the forward model. The location and the mass of the particles are then retrieved by solving an inverse problem, in which the differences between the measurements and the signals generated by the forward model are minimized. This is typically performed by means of appropriate regularization techniques. The MRX forward model and the solution of the inverse problem (imaging) are described in more detail in section 2.1.
An MRXI setup can be broken down into three fundamental constituents: (i) the region of interest (ROI) which is the volume that is imaged and (presumably) contains the MNP distribution, (ii) several excitation coils that produce electromagnetic fields along which the particles are aligned and (iii) a sensor arrangement that records the relaxation signals of the MNPs after the excitation coils have been deactivated. Apart from certain constructional constraints, the sensors and the coils can be arbitrarily positioned around the ROI. To date, a definitive strategy on how to design an MRXI setup in order to obtain the most accurate reconstruction of MNP distributions inside the ROI has not been established. There have been several approaches improving various figures of merit that are designed to predict the reconstruction accuracy based on the assessment of the mathematical model equations of an MRXI system matrix. For instance, these attempts included the improvement of the spatial sensitivity inside the ROI (Coene et al 2012, the maximization of the theoretical information content of the measurements (Coene et al 2015) (see also Krause et al (2008) for a non-MRXI application thereof), the minimization of the condition number of the system matrix (Baumgarten et al 2013, Van Durme et al 2018, Schier et al 2019 as well as-although previously not applied to MRXI improvement but a promising candidate nonetheless-the minimization of the frame potential (Ranieri et al 2014). All of these figures of merit that can assess MRXI systems are described in section 2.2 in more detail.
Naturally, it is desirable to design an MRXI system that achieves the most accurate reconstructions for a given set of constraints (e.g. number of coils, sensors and magnetic field configurations, positions and orientations of coils and sensors, etc), independently of the actual MNP distribution. Furthermore, for biomedical applications, especially for diagnostic and therapeutic purposes, it is crucial to minimize the dataacquisition time of MRXI, which is directly proportional to the number of different applied magnetic field configurations. Each of the previously described approaches aimed to improve the solution of the inverse problem and resulted in enhanced reconstruction qualities for their respective applications. However, it still remains unclear which approach yields the highest benefit for the enhancement of MRXI reconstructions due to the incomparability of the studies' individual results.
Therefore, the aim of this simulation study is to objectively compare the reconstruction accuracies originating from MRXI setups optimized with respect to each of the aforementioned figures of merit. For this purpose, we formulate appropriate cost functions for the figures of merit. The optimizations of these cost functions are then performed using an adapted tabu search algorithm starting from identical MRXI setups and constraints and the reconstructions using the optimized setups are evaluated based on a number of different MNP phantoms. The extensive amount of synthetic data unprecedented in previous MRXI studies that is going to be generated throughout our simulations will for the first time allow a fair assessment of the coherence between the respective figures of merit and the expected reconstruction accuracies. On the one hand, this will lead to an understanding about how the different figures of merit influence the characteristics of the MRXI systems and how reliably their optimization improves the resulting reconstruction accuracy. More importantly, this objective comparison shall also reveal the most effective figure of merit which will play an essential role in the designing process of future MRXI setups.
The optimization is performed by varying the excitation coil positioning and orientations since it already has been shown that MRXI reconstruction qualities heavily depend on the coil positioning (Hoo et al 2014). Also, the positioning of sensors is less flexible compared to coils, for instance due to additional constructional constraints such as a mandatory cooling unit for some types of sensors (Schnabel et al 2004). However, the findings presented in this regard are not limited to the coil configurations but are also applicable to other design parameters of the imaging modality such as sensor positioning and magnetic field configurations. Thus, the results obtained from this study will be valuable for the entire designing process of future MRXI setups.

Methods
2.1. Magnetorelaxometry imaging 2.1.1. Forward model A brief description with the essential modelling steps of the MRXI system matrix L is given below and a graphical representation of the essential model parameters is shown in figure 1. For a more detailed understanding about the modelling process the reader is referred to Crevecoeur et al (2012), Liebl et al (2014).
The ROI containing the MNP distribution is subdivided into N v equally sized cubical or cuboidal volume elements (voxels). A voxel v that contains MNPs with a total iron mass of c v gives rise to a net magnetic moment after magnetization with a previously applied external magnetic field. The magnetic moment is modeled as a point-like source in center Î  r v 3 of the vth voxel, where Î  h v 3 denotes the magnetic field vector at r v . Here, χ denotes the mass magnetic susceptibility of the MNPs and k is a relaxation constant that allows for a stationary formulation of the MRX forward model and which depends on properties of the particle ensemble such as its size distribution and on the properties of the magnetic field. Both constants can be gathered from reference measurements and are therefore known here (see Liebl et al (2014Liebl et al ( , 2015 for more details). The resulting magnetic field emanated by m v can be recorded by N s highly magnetosensitive sensors which are placed at the locations Î  r s 3 . Their respective sensitive axes are denoted by the unit vectors Î  n s 3 . The ultimately measured amplitude of the magnetic flux density B s in sensor s is then given by the superposition of contributions over all N v voxels where μ 0 = 4π · 10 −7 H m −1 is the magnetic permeability of free space and r s,v = r s − r v is the vector between voxel v and sensor s. The operator P · P 2 denotes the Euclidean vector norm. Note that the iron masses c v inside the voxels are the only unknowns of the equation which can be expressed as the vector Î  c N v . The remaining, known parameters in (2) assemble the system matrix Î L a N N s v that encodes the impact of all N v voxels on the N s sensors for a specific magnetic field configuration, referred to as activation sequence a. This enables the formulation of the linear forward model for MRXI representing the recorded relaxation amplitudes of all N s sensors for the activation sequence a. Typically, measurements for MRXI consistd of N a different activation sequences (e.g. by sequentially energizing N a different coils) to increase the stability of the inverse problem and to enhance the reconstruction accuracy of c as a consequence. The multitude of relaxation amplitudes and forward models are concatenated below each other, such that respectively. The complete system of linear equations for the MRXI forward model is then described by Note that, for a fixed voxel and sensor grid, the individual system matrices L L ... N 1 a solely differ from each other by the magnetic field vectors h v,a applied throughout the N a various activation sequences. We sequentially activate each individual coil with unit current in our simulations such that N a equals the number of excitation coils employed in the respective MRXI setups. Consequently, the positions of the electromagnetic coils around the ROI have a major influence on the structure of the system matrix L and thus on the stability and accuracy of the inverse problem's solution.

Solution of the inverse problem
The actual imaging process of MRXI consists of solving the inverse problem of the MRX forward model. For this task, an appropriate regularization technique is mandatory to find a suitable approximation c to the optimization problem with  c ( ) representing the regularization term which is chosen in accordance with the problem at hand. In this study we select two different regularization methods to reconstruct the tested MNP phantoms.
The first method is a classical iterated Tikhonov regularization (Engl et al 1996) with additional nonnegativity constraint (as negative iron masses in voxels are physically not sensible) and is given by where α denotes the regularization parameter. This reconstruction method is giving preference to smooth solutions for c by penalizing its ℓ 2 -norm. The other regularization we apply is the sparsity enforcing iterative shrinkage-thresholding algorithm (ISTA) (Beck and Teboulle 2009), again with added nonnegativity constraint. This method promotes solutions with few nonzero entries in c and is formulated by the optimization problem ( )which are designed to assess the coil subset under consideration regarding its ability to accurately reconstruct the true MNP distribution are introduced in the following (section 2.2.2-2.2.5). Subsequently, optimization problems for the figures of merit are formulated in order to select coil subsets * which produce system matrices  L * ( )that are optimal with respect to the individual figures of merit.

Optimization of spatial sensitivity variation
The spatial sensitivity Î  s N v is measuring the theoretical, overall impact of a voxel on all measurements and has already been used as parameter to increase the performance of MRXI applications . The sensitivity s v for the vth voxel is calculated as the column-sum of absolute values of the system matrix where l i,v is the element in the ith row and the vth column of  L * ( ). It has been shown that a uniform spatial sensitivity distribution acts beneficial on the reconstruction accuracy , Coene et al 2015. Therefore, we choose an objective function that quantifies the variation of s and seek to minimize it. The coefficient of variation which is the ratio between standard deviation STD and mean μ, measures the dispersion of s. It is employed here as figure of merit to pick a coil subset optimal with respect to the spatial sensitivity distribution. Thus, the corresponding optimization problem is formulated as . 10 c * * * * * ( ( ( ))) | | ( )

Optimization of mutual information
The mutual information MI is another parameter that has been used before to improve the system matrix with respect to solving the inverse problem of MRXI (Coene et al 2015). The mutual information of two random variables X and Y is defined as the difference between the entropy of one variable H(X) and the conditional entropy of both variables H(X|Y) (Cover and Thomas 2012) MI is commonly used in information theory and can be interpreted as the amount of information X and Y have in common. In our case, the individual random variables consist of the system matrices Î L i N N s v with Î  i that arise from solitary activation of coil i with unit current. The normalized mutual information scales the MI to a range between 0 and 1. This is especially useful when evaluating more than two random variables, which is the case when multiple different system matrices L i constitute the final system matrix  L * ( ). The goal is to select a coil subset *, where each coil introduces as much new information into the system equations as possible. This leads to an optimization problem where we want to minimize the sum of NMI over all possible combinations of L i in the subset * such that

Optimization of condition number
The decrease of singular values s s s = L L L ...
with l i and l j denoting the ith and jth rows of  L * ( ). In frame theory, FP is interpreted as the potential energy between the constituent vectors of a frame (Casazza et al 2006). It is large if the angles between the vectors are small (or close to integer multiples of π). Conversely, the FP is low if the vectors are optimally spread on  N v such that the angles are as large and as uniform as possible which implies an ideal coverage of the search space. According to frame theory, it has been shown that a system matrix with minimal FP theoretically yields the lowest mean square error between the least squares solution c and the true image c (Benedetto and Fickus 2003).
Therefore, an objective function seeking a coil subset with minimal FP is formulated:

Optimization algorithm
The number of different coil combinations N c * out of a larger coil configuration with N c coils is given by Even for a relatively small coil configuration with N c = 100 and a small subset of = N 10 c * coils, this leads to an order of approximately 10 13 combinations. Clearly, a direct computation of all possible choices is not feasible for larger coil configurations. We employ a global tabu search optimization algorithm to determine subsets * that result in minimal values  f L * ( ( )) of the figures of merit introduced in section 2.2. Presented by Glover (1990), Glover and Laguna (1998), the algorithm was already employed for the sensor selection problem in Lau et al (2008), Eichardt (2012) and is adapted to the coil selection problem in the following.
The proposed tabu search algorithm is summarized by algorithm 1. The input parameters consist of the system matrix arising from individual coil activations of all possible coils  L( ), the number of coils in the entire coil configuration N c , the number of coils to select for the subset N c *, the minimum and maximum number of neighbours to test < n n min max , as well as the list  . This list contains the numbers of coils of the current best solution * that are to be randomly exchanged with unselected coils from  * ⧹ during each for-loop iteration of the algorithm in a descending order with " The basic idea of algorithm 1 is to start with a globally oriented optimization approach, where a large number of coils (defined by  1 [ ]) in a randomly initialized subset * are exchanged with randomly selected unused coils from the set  * ⧹ . This exchange is controlled by the function ) which then serves as new starting point for the next coil exchange. The number of new candidate solutions n (i.e. coil exchanges to perform) that shall be tested in the first iteration of the for-loop i = 1 is given by = n n min . Each time a new candidate solution is tested, it is added to the tabu list , thereby avoiding evaluating identical subsets repeatedly. Throughout the iterations of the for-loop, the number of coils to replace  i [ ] is progressively decreased while the number of new candidate solutions n is progressively increased (controlled by the function ¬ n getNumNeighbours i n n , , min max ( )) in order to transition from a globally oriented to a locally oriented optimization approach. The getNumNeighbours function interpolates linearly between n min and n max with  | | steps and returns the ith interpolated point rounded to an integer value. Finally, in the last iteration of the forloop =  i | |, the largest number of = n n max new candidate solutions are evaluated.

Setups and phantoms
Two different MRXI setups for evaluating the optimization strategies are employed to determine if the preferential figures of merit for coil positioning are consistent over alterations in the setup design such as number and size of voxels, shape and resolution of the ROI, as well as varying sensor distance. The sensor grid remains unchanged for both setups. The sensor positions and orientations employed in our simulations are selected according to the real 304 SQUID (Superconducting Quantum Interference Device) magnetometer grid used at PTB Berlin (Schnabel et al 2004) since it is commonly applied in numerous MRXI studies and experiments. The 304 sensors have been arranged in four horizontal planes and are oriented in different directions (see blue arrows in figures 2(a) and (b)).
The first setup (setup 1) employed here consists of a cuboidal ROI with the dimensions 12 × 12 × 6 cm 3 which has already been used in other MRXI studies , Liebl et al 2015 (see figures 2(a) and (c)). The ROI is located centrally below the sensor system with a distance of 6 cm between the lowest sensor plane and the upper boundary of the ROI. It has been tesselated into a 10 × 10 × 5 voxel grid to enable the formulation of the forward model (4). A total of 210 electromagnetic coils are placed on four sides of the the ROI in negative and positive x and z directions, henceforth referred to as x − , x + , z − and z + . The other two sides need to be accessible for the experiments such that no coils can be placed there. Regular coil configurations are arranged in six planes (p 1 -p 6 ) in each direction, respectively, except for z − , where two additional planes holding large coils have been placed (p 7 and p 8 , see figure 2(c)). In general, flat spiral coils with a wire distance of 0.127 mm have been employed in this study. Only on plane p 5 , three additional, larger circular coil loops have been placed around the spiral coils to allow for an increase in coil radius without having to increase the distance to the ROI (typically, closer coils provide a better spatial encoding of the ROI). The coil radii range from 1 cm on p 1 up to 10 cm on p 8 and exhibit interplanar distances of 0.5 cm between the planes p 1 -p 5 and 1 cm between the planes p 5 -p 8 (the distance between the more distant planes has been increased to avoid overlap of the coils). The computed magnetic field amplitudes in the center of the ROI range between 0.2 and 157 A m −1 .
Setup 2 has a tomography styled design layout with a cylindrical ROI that measures 20 cm along its central axis in y direction and has a radius of 5 cm (see figure 2(b)). The ROI lies centrally beneath the sensor grid and the vertical (z) distance from the border of the ROI to the lowest sensor plane is 5 cm. Initially, it has been discretized into 10 × 10 × 10 voxels on the basis of a tight cuboidal volume around the cylindrical ROI. However, all volume center points that lie outside of the ROI are neglected in the simulations, leaving a total of 800 voxels in an approximately cylindrical shape. In sum 519 electromagnetic coils are arranged in six tubular regular configurations with increasing radii (t 1 -t 6 ) around the ROI (see figure 2(d)). The narrowest tube (t 1 ) exhibits a radius of 5.5 cm which increments by 0.5 cm for each of the grids t 2 -t 4 . The radii of the grids t 5 and t 6 measure 7.75 cm and 9.5 cm, respectively. The spiral wire coil's center points coincide with the six tubes and are oriented towards the central axis of the cylindrical ROI. The smallest coils, which are closest to the ROI (t 1 ), exhibit coil radii of 0.75 cm and the largest coils at t 6 measure 8 cm in radius. The computed magnetic field amplitudes in the center of the ROI range between 0.006 and 1.7 A m −1 .
We use three distinct MNP distributions (i.e. phantoms) to evaluate the reconstruction accuracies of the different coil selection strategies over varying shapes of the MNP ensembles. The phantoms are shown in figure 3. Each phantom exhibits a different property that can be an obstacle with respect to accurate MRXI reconstruction results. Phantom P consists of a sparse MNP distribution in form of the letter 'P' in the center of the ROI. Here, the difficulty lies within the reconstruction of the cavity in the center of the phantom. The phantom PTB is composed of its eponymous letters as well, where 'T' is located in the central plane and 'P' as well as 'B' are located in the outermost planes of the ROI. The challenge of this phantom is the reconstruction of three different MNP accumulations in spatially separated domains. The planes phantom consists of five large cuboidal MNP clusters that are arranged edge to edge in an alternating manner as shown in the right panels of figure 3. This phantom is typically complicated to reconstruct due the strong spatial variations of MNP concentrations in the voxels. The particle characteristics are taken from the experiments conducted by Liebl et al (2015) with an iron concentration of 3.7 mg ml −1 and a product of relaxation constant and mass magnetic susceptibility of approximately k · χ = 6 · 10 −4 m 3 kg −1 (known from reference measurements).
2.5. Coil subset evaluation 2.5.1. Assessing the figures of merit The tabu search (see section 2.3) has been carried out for multiple subsets with different numbers of coils which empirically yielded a reasonable trade-off between computational effort and algorithm convergence.
The optimization is repeated 10 times for each element in  c , producing 10 (potentially) different subsets for every N c * tested in both setups which resulted in a total of 3520 tested coil subsets over both setups, all figures of merit and all repetitions of the tabu search. Measurement data b is simulated based on every single subset for all three phantoms, where each coil is individually driven by unit current. Prior to the reconstruction of the MNP distributions, white Gaussian noise with SNR = 20 dB is added to every measurement b employed in this study. This noise level is commonly observed in MRXI studies Haueisen 2010, Coene et al 2014). Every measurement b is corrupted by 10 different, newly generated white Gaussian noise patterns, thereby realizing 10 slightly varying and noisy measurements. The reconstructions have been performed based on these 10 different measurements for all subsets and each of the three phantoms. The accuracies of the reconstructions are assessed via the Pearson correlation coefficient CC (Pearson 1896, Dunn et al 1987, where the reconstructed and the ground truth MNP distributions are compared and the similarity between the two images is quantified in a range from −1 (perfect anticorrelation) over 0 (no similarity) to 1 (perfect reconstruction). The reconstructions of every single measurement are repeated with a variety of different regularization parameters α (see (6) and (7)). α is chosen in the range 10 −6 α 10 −1 in 11 logarithmically equidistant steps for the Tikhonov regularization and in the range 10 −4 α 1 in 9 logarithmically equidistant steps for the ISTA regularization. These boundaries were empirically determined such that every reconstruction result progressively deteriorated with any further reduction/increase of the smallest/largest regularization parameter values. For every measurement, the reconstructions with the maximum CC values achievable through this variation of the α values are stored for the subsequent evaluation of the reconstruction performances. The mean values and the standard deviations of the joint 100 maximum CC values from the 10 measurements corrupted by different white Gaussian noise patterns for each of the 10 (potentially) different subsets have been calculated for every coil count N c *. This is done separately for each of the four figures of merit, all three phantoms and both setups. Hence, a total of 105 600 reconstructions has been ultimately considered for the evaluation of the coil subset selection strategies.

Numerical results
We compare the coil subsets that are selected with respect to the four different figures of merit presented in section 2.2 regarding their ability to accurately reconstruct different phantoms. Figure 4 depicts exemplary specimens of coil subsets for both setups that were selected using the tabu search algorithm (see section 2.3) with respect to all figures of merit and a coil count of = N 20 c * , respectively. Generally, the subsets produced by the four optimization problems differ strongly from each other by visual inspection. On the other hand, subsets selected with respect to a specific figure of merit show similar positioning tendencies. Hence, the subsets of each optimization problem exhibit individual characteristics which are explained in the following.
The optimization with respect to the coefficient of variation of the spatial sensitivity CV(s) preferably selects larger coils as they provide a more uniformly distributed sensitivity throughout the ROI. Furthermore, smaller coils at the corners and edges of the ROI are also preferred for the selection to compensate the lower sensitivity values of the larger coils in these areas. Additionally, the optimization chooses the coils in z − -direction rather than on the opposite side of the ROI to resolve the lower sensitivity because of the increased distance between the voxels and the sensors. The tabu search employing the normalized mutual information NMI as objective function seeks to maximize the theoretical information content of the system matrix by choosing a large variety of different coil diameters and locations. When optimizing the condition number κ, smaller coils are preferred over larger coils, especially if they are located on top of the ROIs (z + ). Similar to the optimization of CV(s), the focus lies here on selecting coils in the lower regions (z − ) as well. However, in contrast to the subsets minimizing CV(s), the chosen coils for optimization of κ are mostly positioned centrally beside the faces of the ROI, specifically for setup 1. The subsets yielded by the optimization of the frame potential FP almost exclusively employ the smallest coils available for selection. Additionally, the coils are always well distributed over the entire ROI. Figure 5(a) shows the CC mean values and standard deviations that were reached with coil subsets from setup 1 using coil numbers   N 1 64 c * as described in section 2.5.1. The reconstructions are performed with Tikhonov regularization (top row) and ISTA regularization (bottom row) for the phantoms P (left column), PTB (middle column), and Planes (right column). Reconstructions of subsets with coil counts  N 5 c * (gray areas in figure 5) are very volatile such that no meaningful comparison can be made and are therefore neglected in the present evaluation. For larger coil counts, coil subsets selected with respect to minimal condition number κ are among the highest CC values in almost all cases. The subsets selected according to FP also produce good reconstruction results in most scenarios, however, they are not able to surpass the κ-subsets, save for some rare exceptions. NMI is the third best figure of merit for setup 1, already with considerable deterioration in reconstruction accuracy compared to κ and FP. Additionally, subsets optimized with respect to NMI yield inconsistent reconstruction accuracies, indicated by the large standard deviations of the CC values. The least efficient objective function is the sensitivity variation CV(s) which yields far worse reconstruction results than any other objective function in all cases for setup 1.
The reconstruction results with optimized subsets of setup 2 are shown in figure 5(b) and have been performed under the same conditions as for setup 1. The best performance indicator out of the four figures of merit is not as distinctly classifiable as for setup 1. However, coil subsets selected with respect to FP yield worse reconstruction results compared to the other three types of subsets in most cases. On average, the other three figures of merit, namely κ, NMI as well as CV(s), produce similar CC values. While the subsets selected via κ and NMI yield almost identical reconstruction accuracies in all cases, the CC values produced by CV(s)-subsets are phantom dependent. Finally, both κand FP-subsets produce inconsistent imaging results in most reconstruction scenarios of setup 2 which is indicated by their enlarged CC value standard deviations.
Additionally, the quantification performances of the individual subsets were assessed using the root-meansquare deviations RMSD between the reconstructions and the true MNP phantoms. It is worth mentioning that the α values that produce the largest CC values also yield the lowest RMSD in almost all cases. Exemplary quantification results from the different setups, phantoms and regularization techniques over all employed N c * are shown in figure 6. In general, the quantification performances show trends very similar to the reconstruction accuracies depicted in figure 5. Thus, the conclusions that can be drawn from the quantification performance assessment are closely related to the ones deduced from the reconstruction accuracy evaluation.
On average, the subsets chosen according to the condition number κ yield the highest CC values as well as the lowest RMSD values when taking all reconstruction results from both setups into account. The κ-subsets converge within approximately 20-30 coils per subset to the largest achievable reconstruction accuracies of all coil configuration selection strategies for setup 1. This coil count is in line with the findings reported in Coene et al (2014Coene et al ( , 2015. Due to the finer voxel grid, the CC values of κ-subsets do not plateau within the tested coil counts of  N 64 c * for setup 2, but are among the best reconstruction results in most cases nonetheless. Almost identical observations are made in the quantification performance evaluation. In setup 1, FP shows a rather good reconstruction performance as well, however, the CC values achieved with this figure of merit in setup 2 are mostly far below the other strategies. Analogously, NMI and CV(s) yield results similar to those of κ in setup 2, but strongly lack reconstruction accuracy in setup 1 in comparison. Like before, the assessment of the quantification performances shows a very similar behavior here. This leads to the conclusion that coil subsets selected with respect to κ are the most reliable to offer an accurate reconstruction of the MNP distributions. Thus, the coil configurations optimized via κ are further analyzed in the following.
The coherence between the tested figures of merit and the achievable CC values is evaluated based on all optimized subsets of setup 1 with = N 20 c * coils. Setup 1 has been used due to the fact that the CC values of the different subsets are spread over a larger range than in setup 2 which simplifies the detection of a relation between the figures of merit and the achievable reconstruction accuracy. = N 20 c * also has been chosen because of the comparably large differences between the CC values at this coil count. The mean values and standard deviations of all figures of merit are calculated over the 10 subsets obtained from each of the four optimization strategies, respectively. Analogously, the average CC values and their standard deviations are calculated over the total of the 100 reconstructions for every optimization strategy. Each row in table 1 shows these values for the subsets optimized with respect to one of the four different figures of merit. A clear inverse relationship between the condition number κ and the resulting CC values is visible in the data presented in the table. However, no distinct correlation between reconstruction accuracy and the other three figures of merit can be determined from this evaluation. The κ-subsets, which produce the highest CC values, are on average only second best in FP, third best in CV(s) and yield the worst NMI values.
A similar result is visible in the evaluation depicted in figure 7. Here, the figures of merit of the three subsets optimized with respect to κ that yield the highest CC values in setup 2 are averaged and compared to the corresponding averages of the three κ-subsets with the worst reconstruction performance. This procedure has been repeated for every coil count N c *. We use setup 2 for this evaluation due to the larger standard deviations of the CC values of κ-subsets that indicate a larger performance gap between the achievable reconstruction accuracies among the subsets than in setup 1. This allows for a clearer identification of a possible relation between the figures of merit and the CC values. On average, the better performing subsets exhibit lower condition numbers than their worse performing counterparts for almost every coil count N c *. Apart from that, the other figures of merit give no conclusive indication whether or not a coil subset is well suited for accurate reconstructions of MNP distributions. Both CV(s) and NMI fluctuate randomly over the examined coil counts N c *. The majority of better performing subsets even show a slightly increased frame potential FP which is in direct contradiction to the original assumption that a lower FP indicates higher reconstruction accuracy.
Finally, the achievable reconstruction accuracy of an exemplary κ-subset with = N 20 c * is compared to the reconstructions obtained from a regular coil configuration using setup 1 and the PTB phantom in figure 8. This is done to demonstrate the improvement of the reconstruction accuracy of coil configurations optimized with respect to the condition number compared to an intuitive coil configuration design. The regular coil configuration has been assembled such that the four available faces of the ROI are well covered by five excitation Figure 6. Exemplary quantification performances of the differently optimized setups, phantoms and regularization techniques, analogous to the evaluation shown in figure 5. The quantification performances are measured by the RMSD between reconstructions and ground truth phantoms. Again, the circles indicate the average RMSD values, the pastel-colored backgrounds depict the respective standard deviations, and the datapoints in the gray areas (  N 5 c * ) are neglected for the comparison of the results. coils, respectively, as depicted in the top left panel of figure 8. The reconstructions have been performed based on sequential coil activations using the Tikhonov regularization and the resulting images with the largest achievable CC values are shown below the corresponding setups in figure 8. It is possible to recover the PTB phantom with an accuracy of CC = 0.81 using the regular configuration. In comparison, the κ-subset yielded a reconstruction with an accuracy of CC = 0.95. Hence, a far more accurate reconstruction of an underlying MNP distribution using the same amount of coils/activation sequences is possible by means of a system matrix optimized with respect to κ compared to an intuitive design of the coil configuration.

Discussion
We assess various strategies for optimal positioning of excitation coils as well as optimal coil radius selection for MRXI setups in this extensive simulation study. Here, the goal is to find coil arrangements for a given ROI, sensor setup and number of coils N c * which enable the most accurate reconstructions of different MNP distributions using a vast amount of synthetic data unprecedented in previous MRXI simulation studies. To this end, two different MRXI setups have been generated, each of which employs a large amount of excitation coils (see figure 2). A global optimization algorithm in the form of a tabu search (see section 2.3) has been applied to determine subsets of coils with a predefined desired coil count N c * with respect to several different figures of merit. In a prior step, these figures of merit have been implemented into corresponding cost functions, thereby enabling their direct optimization. The applied figures of merit evaluate and quantify certain properties of the system matrices resulting from individual, sequential activations of the coils in the subsets driven by unit current. All four of these examined figures of merit shall provide a prediction about whether a system matrix is well suited for the reconstruction of an MNP distribution, or not. This allows for the first time an objective comparison between the different figures of merit regarding their ability to predict MRXI reconstruction accuracy.
The coil configurations optimized in this study with respect to the individual figures of merit are largely in accordance with the findings reported in literature for manually or randomly determined coil positions. Furthermore, the results presented in this study provide new insight as to why certain outcomes were previously observed in literature. For instance, the coils of κ-subsets are mainly concentrated below the ROI and become scarcer the closer they lie to the sensor grid. Additionally, larger coil radii are mainly found below the ROI. This coincides with the optimized MRXI coil currents reported in Schier et al (2019), where the lower coils are driven by much larger coil currents than the upper coils. Also the selection of the coil radii according to κ is reasonable since smaller coils have been reported to generally reconstruct more accurately than larger coils (Hoo et al 2014). The subsets selected via CV(s) show an increased concentration of coils in the lower regions of the ROI as well. The coil current patterns determined by Crevecoeur et al (2012), who also sought to equalize the spatial sensitivity distributions, show the same behaviour-namely, increased driving coil currents of the lower coils in their setup. The coil configurations optimized via NMI are well distributed around the entire ROI and employ many different coil radii. Although there are no comparable results regarding the selected coil radii in literature, the coil positioning presented here is largely in accordance with Coene et al (2015) (which does not display the actual coil configuration but refers to an analogous coil configuration presented in Coene et al (2014)). The results reported in Coene et al (2014) also indicate that a well distributed coil positioning is beneficial for low NMI values. Regarding the use of FP, the benefit of selecting small coil radii according to Hoo et al (2014) was demonstrated. Other studies on this topic focus on the quality of their solutions and not on their optimized variables (Joshi and Boyd 2008, Ranieri et al 2014, Jiang et al 2016. We formulated corresponding cost functions for each of the four figures of merit (see section 2.2) with the goal to optimize MRXI setups and to improve the solution accuracies of the underlying linear inverse problems. This allowed for a fair comparison among the optimization results, thereby enabling the determination of the optimization approach that leads to the largest and most reliable improvement of MRXI reconstruction accuracy. In this study, we selected several coil subsets optimal with respect to each of the four figures of merit with coil counts   N 1 64 c * , respectively, thereby simultaneously optimizing coil positioning and sizes for the first time in MRXI. Two different regularization approaches have been tested to reconstruct three phantoms with each of the subsets. The corresponding reconstruction accuracies (measured via CC) are displayed in figures 5(a) and (b). Additionally, the quantification accuracies were assessed by calculating the respective RMSD values (see figure 6 for exemplary quantification results). These results show that, on average, the subsets selected via the condition number k  L * ( ( )) enable the most accurate reconstructions of the phantoms in the majority of the cases, or are at least close to the best performing subsets. This is evident in both tested MRXI setups, which indicates that the optimization approach reliably promotes good imaging qualities independent of the employed shape and resolution of the ROI (i.e. number, size, and location of the voxels). However, one also needs to keep in mind that these coil configurations are optimized for specific voxel grids and that the imaging performance might vary when changing the shape or the resolution of the ROI (although the resolution has a much smaller impact on the optimization result). An additional optimization procedure could be necessary, depending on the alterations that are required.
The subsets optimized with respect to the other three figures of merit yield inconsistent reconstruction performances across the two different setups. The evaluations presented in table 1 and in figure 7 support this conclusion. Here, the only reliable relation between the figures of merit and the reconstruction accuracy is that lower condition numbers consistently yield higher CC values. Thus, the minimization of k L ( ) can be considered as the most reliable and efficient approach for optimizing an MRXI setup if no a priori information about the MNP distribution inside the ROI is available.
Nonetheless, the reconstruction results can be enhanced further if the approximate location of the MNPs is known beforehand. As an example, the subsets optimized with respect to CV(s) yield the most accurate reconstructions in setup 2 for the PTB phantom. This is the case because, as mentioned in the results section, the tabu search optimizing CV(s) preferentially selects coils at the edges of the ROI. Since a large portion of MNPs in the PTB phantom is concentrated in regions close to the edges, these selected coils provide a considerable information gain in those areas. This effect was also observed in Coene et al (2012).
Based on their large standard deviations the selection strategies show inconsistent reconstruction results in some cases, such as NMI in setup 1 and κ as well as FP in setup 2. In case of NMI this can be attributed to the figure of merit itself. The subsets selected via NMI are always constituted of both very large as well as smaller coil radii in order to maximize the theoretical information content in the system matrix. Due to the coarser voxel grid and smaller coil to voxel distances in setup 1 more nonredundant information can be introduced into the system matrix by using larger coil radii than in setup 2. However, in some cases, these larger coils are placed on top of the cuboidal ROI (z + -direction) which introduces a large concentration of spatial sensitivity in the upper regions due to the vicinity to the sensor grid. This destabilizes the inverse problem and causes worse reconstruction results than in subsets where no large coils are positioned on top of the ROI.
On the other hand, the variations in the reconstruction accuracies of κ and FP in setup 2 stem from the vast increase of possible coil combinations due to the larger number of coils to select in combination with the multitude of local minima in these figures of merit in which the optimization algorithm can get stuck. In case of κ, this is visible in figure 7 where the coil subsets that were able to reach lower condition numbers also produce more accurate reconstructions. An increase of candidate solutions n min and n max in the tabu search would result in more consistent coil subsets and thus would decrease the variations in the reconstruction accuracies.
There are also some limitations to the presented approaches and results that have to be addressed. For instance, the numbers of coils in the complete coil configurations of the tested setups have been deliberately chosen relatively small to achieve a good trade-off between the necessary computation time and the quality of the solutions. Naturally, denser coil arrangements employing a multitude of overlapping coils with different coil radii and coil orientations per coil position might lead to optimization results with even lower κ values. However, apart from the overlap constraint that would have to be implemented in the tabu search algorithm, this would vastly increase the required computation time of the presented approach (and also any other global optimization technique) to determine a good minimum of the objective function and would also introduce many more local minima to the cost functions which would render the problem intractable. Therefore, as mentioned before, a continuous local optimization method subsequent to a rough global optimization is much more promising in this regard. Nonetheless, the computed subsets already result in low condition numbers, are able to produce accurate reconstructions of the underlying MNP distributions and provide a good initialization for a local optimization approach.
Obviously, one also has to consider the practicality of the optimization results of the presented approaches. Whereas small coils are certainly beneficial for the minimization of κ due to their highly divergent magnetic fields, they produce very small field strengths (within a certain range of available amperage) which typically mainly sensitize superficial regions of the ROI and eventually produce rather small measurement signals as a result. These small coils are selected in our optimization since κ is based on the mathematical structure of  L * ( ) and does not consider the noise level on the measurements that is to be expected. Theoretically, these coils might be the best choice for completely noiseless measurements, however, in practice sufficiently large sensor signals must be recorded to overcome the ambient noise. In our simulations the effect on the reconstructions is attenuated due to the fact that the noise corruption on the signals is calculated relative to the recorded measurement amplitudes with SNR = 20 dB (i.e. smaller signals are corrupted by weaker noise). However, in practice one is typically confronted with a specific noise profile inherent to the measurement environment. Thus, in further investigations the initial complete coil configurations should be constructed with these considerations in mind (e.g. lower coil radii limits or minimal signal strengths). Also, the expected measurement noise should be integrated already as additional criterion that influences the coil selection process in future optimization approaches. The employed noise model for the simulated measurement data is not entirely in accordance with the noise in practical MRXI measurements. Although the average noise level of SNR = 20 dB coincides with realistic MRXI measurements using similar MNP concentrations and setup dimensions Haueisen 2010, Coene et al 2014), the assumption of a normal distribution for the noise is an approximation. For instance, the SQUID sensor system of PTB Berlin shows several different noise distributions that depend on the orientations of the sensors (Coene et al 2014). Furthermore, deviations in the locations and orientations of the coils and the sensors from the estimated positions used for computing the forward model introduce relatively large model inaccuracies in the system equations. Additional model errors can be related to variations in the magnetic properties of the particles. Among others, this might include varying relaxation characteristics due to MNP aggregation/agglomeration, possible magnetic interactions between the particles, and the slightly nonlinear magnetic field dependence of the relaxation constant. All of these inaccuracies in the applied system equations amount to noise conditions that are hardly predictable and thus have to be approximated. This introduces small uncertainties in the credibility of the reconstruction results presented in this study. However, it is not expected that slight variations of the employed noise distributions would considerably alter this study's simulation outcomes. Practical measurements or a more realistic noise model for the simulations will be necessary to fully validate the presented reconstruction results.
It is also noteworthy that the presented optimizations were performed based on the SQUID sensor system of PTB Berlin that is situated inside the Berlin Magnetically Shielded Room 2 (BMSR-2), which is currently the globally best performing walk-in sensor-shielding combination in terms of sensitivity, shielding factor, as well as number of sensors. Naturally, when employing a different sensor system, considerable changes in the sensor positions and orientations would necessitate a repetition of the optimization procedure that, again, would result in favorable system matrices in a purely mathematical sense. However, also worse noise conditions are expected if an optimized MRXI setup is operated in an environment other than the BMSR-2. This would most notably impact the utility of small coils (due to their comparably small magnetic fields that lead to weak relaxation signals) which further elucidates the requirement of implementing a realistic noise model-including the previously discussed influences of positioning errors and deviations in the magnetic particle properties-into the optimization procedure that actively influences the coil selections based on expected SNR values. The next step for the optimization of MRXI setups is going to be the simultaneous optimization of multiple design parameters such as the coil positioning, radii and orientations, the sensor positioning and orientations as well as the driving coil currents with respect to κ in order to obtain a better optimization result. Especially the optimization of the sensor locations and orientations recently became increasingly attractive due to the emerge of optically pumped magnetometers which are more flexibly positionable than the spatially constrained SQUIDs considered in the present study (Jaufenthaler et al 2020). Since all of these design parameters are not independent from each other, a consecutive optimization of each variable would not be sufficient to determine a good minimum of the objective function.

Conclusion
This simulation study shows that there are considerable differences in reconstruction quality depending on coil positioning and coil radii, clearly visible in the reconstructions depicted in figure 8. Throughout numerous MRXI reconstructions, it has been demonstrated that the minimization of κ is the most promising approach to enhance the accuracy of the solution of the inverse problem, outperforming the other figures of merit (CV(s), NMI, FP) that were applied in literature. Thus, future steps toward the optimal design of an MRXI setup should positively include considerations regarding the effects of the design parameters on the condition number of the system matrix and as such pave the way towards obtaining clinically relevant reconstructions.