Quantitative predictions in small-animal X-ray fluorescence tomography

: X-ray fluorescence (XRF) tomography from nanoparticles (NPs) shows promise for high-spatial-resolution molecular imaging in small-animals. Quantitative reconstruction algorithms aim to reconstruct the true distribution of NPs inside the small-animal, but so far there has been no feasible way to predict signal levels or evaluate the accuracy of reconstructions in realistic scenarios. Here we present a GPU-based computational model for small-animal XRF tomography. The unique combination of a highly accelerated Monte Carlo tool combined with an accurate small-animal phantom allows unprecedented realistic full-body simulations. We use this model to simulate our experimental system to evaluate the quantitative performance and accuracy of our reconstruction algorithms on large-scale organs as well as mm-sized tumors. Furthermore, we predict the detection limits for sub-mm tumors at realistic NP concentrations. The computational model will be a valuable tool for optimizing next-generation experimental arrangements and reconstruction algorithms.


Introduction
X-ray fluorescence (XRF) tomography using nanoparticles (NPs) is an emerging bioimaging modality with potential for improved molecular imaging. Recent experiments with spectrally matched Mo NPs show that such imaging can be performed with sub-200-µm spatial resolution in mice [1]. However, several factors limit the observability and quantitative accuracy in reconstructed images. In the present paper we describe an efficient computational tool that accurately models small-animal XRF tomography in order to provide quantitative and realistic predictions of, e.g., small-feature observability or reconstruction accuracy in 3D.
XRF tomography has long been used at synchrotron-radiation sources for sensitive elemental 3D nano-and micro-imaging [2]. During the last decade the method has been extended to laboratory sources combined with NPs with the goal to provide molecular or functional bioimaging with new contrast and improved resolution. Data acquisition can be based either on a scanning pencil beam geometry [1,[3][4][5][6], cone-beam geometry with a collimated detector [7][8][9] or parallel approaches with single or multiple pinholes on the detector side [10,11]. The pencil-beam approach typically yields higher spatial resolution while the parallel methods have potential for faster data acquisition. Experimentally, these methods have primarily been tested on phantoms. Recently, however, the pencil-beam approach has been demonstrated in mice models with gold NPs [9] with low resolution and molybdenum NPs with high resolution [1] in laboratory arrangements, and earlier at a synchrotron studying cerebral iodine perfusion [12].
Simulations play a key role in predicting the performance of XRF tomography arrangements [7,[13][14][15][16]. While observability largely relies on experimental factors such as radiation dose, NP concentration and spatial resolution, the quantitative accuracy of the tomographic reconstructions relies heavily on a correct physics model when reconstructing the data. Evaluation of the latter has so far been limited to either overly simplified phantoms with a priori information about the composition, or complex experimental samples with no prior information, making quantitative estimates of reconstruction accuracy difficult.
Such simulations are preferably based on Monte Carlo (MC) approaches. Many MC codes exist today and some have been used to simulate XRF tomography (PENELOPE [13], Geant4 [14], FLUKA [15], EGSnrc [16]). While these general-purpose tools are well-validated, they all suffer from computational inefficiency. They are usually designed to run on single CPUs with little to no support for parallelization, leading to very long simulation times. Even on a high-end computer, single slice XRF tomographic simulations can easily take up to a week [16]. This has forced the use of large computational CPU clusters in order to make simulations feasible time-wise [14,15]. The lack of computational efficiency has resulted in previous studies being limited to single slice, i.e., 2D imaging, typically of oversimplified phantoms.
Graphics processing units (GPUs) are increasingly used for scientific computation and a number of GPU-accelerated Monte-Carlo tools for x-ray transport have emerged [17][18][19]. The architecture of the GPU allows simulating photons in parallel, leading to substantial accelerations in computational speed. Typically, acceleration is achieved through clever utilization of the different GPU memory levels in combination with a stripped physics model. The latter decreases the number of possible paths a simulated photon can take, making the acceleration more efficient. However, this also means that the codes become application specific, since the stripped physics model is not suitable for all applications. As an example, the accelerated MC-GPU code [17] simulating cone-beam computed tomography (CBCT), completely neglects secondary emission of electrons and XRF photons.
In the present paper we describe a GPU-accelerated Monte-Carlo simulation tool which is up to 3 orders of magnitude faster than tools previously used for simulating XRF tomography. The high speed allows 3D simulations of, e.g., detailed small-animal phantoms providing accurate quantitative predictions of small-animal XRF tomography. Furthermore, we can quantitatively evaluate the physics model in the reconstruction algorithms in realistic smallanimal imaging scenarios. We demonstrate the tool on our laboratory pencil-beam arrangement, both for large-scale (lungs, liver etc.) and small-tumor quantitative imaging, as well as predictions on sub-mm tumor detectability. In summary, the unique combination of a highly accelerated MC tool combined with accurate small-animal phantoms allows unprecedented realistic full-body tomographic simulations that can be performed on a single computer.

Computational model
We have developed a GPU-based XRF tomography simulation tool, "XRF-GPU". To our knowledge, it is the fastest MC tool for simulating XRF tomography. It is capable of simulating practically any voxelized object and data acquisition geometry, with fully customizable source and detector configurations. An early version [20] was used to study different anti-scatter grid solutions for XRF tomography [21].
The code is based on the open source MC-GPU [17]. In addition to the features included in the original code, we have implemented the physics necessary to simulate XRF. Each time a photon is absorbed in a certain region of the simulated object, a function is called to see if the deposited energy is high enough to ionize inner-shell electrons in the region. If true, the incident photon is transformed into an XRF photon and emitted isotropically. Since the material of the absorbing region is typically composed of several different elements, first the absorbing element is determined using probabilities from the a priori known weighted elemental composition of the material. Then, the energy of the XRF photon is chosen according to tabulated probabilities for α K and β K emission for that specific element [22].
The speed of XRF-GPU enables 3D simulations of high-resolution voxelized objects, which allows us to combine it with available small-animal digital phantoms to create a realistic preclinical computational model. We create a voxelized mouse phantom from the open source D organs segme organ tissues artificially "in shows the mo realistic imag used in the tomography ( We perfor resulting in a with the real tumors with g in 8 hours usi Extrapolating that with XR projection cou a single deskt

Reconstruction algorithms and self-absorption correction
To achieve quantitative reconstructions, we correct for self-absorption by using the transmission sinograms recorded simultaneously as the XRF sinograms. We reconstruct the transmission sinograms using filtered back-projection (FBP) to achieve a set of CT images. These can also be seen as absorption maps recorded at 24.1 keV, and are therefore used to create corresponding absorption maps at the 17.4 keV Mo α K XRF. Since we image smallanimals, we do this by making the dual material approximation of soft-tissue and bone and then use the corresponding tabulated attenuation ratios for each respective region [30]. This absorption map is then used to correct for the XRF self-absorption [1].
We capture the physics of the XRF imaging process in a quantitative forward model. It consists of the XRF self-absorption map together with further a priori information about the imaging system such as incident flux, detection solid angle, detector efficiency etc. XRF sinograms are pre-processed in order to estimate and remove background contribution from Compton scattering. Estimation is done by linear interpolation of the photon counts on each side of the XRF peak. The pre-processed XRF sinograms are then fed, together with the forward model, into our simultaneous iterative reconstruction technique (SIRT) based reconstruction algorithm to achieve quantitative 3D reconstructions of the NP distribution inside the small-animal.

Quantitative large-organ simulations
XRF tomography experiments show significant accumulation of NPs in major organs such as liver, lungs and kidneys [1,9]. While this is typically unwanted in the case of tumor localization, it is important to quantify the amount of NPs that has been accumulated in these organs. This poses an interesting reconstruction challenge. In the large volumes of the major organs, the region closer to the core of the animal will experience much more severe selfabsorption than regions close to the skin. Thus, a homogenous distribution in a large organ may appear inhomogeneous in the case of an improper self-absorption correction. The XRF-GPU computational model allows us to assess how well our algorithms can reconstruct the homogeneity at different concentrations.
Using our mouse phantom, we artificially "inject" the lungs with homogeneous Mo NP concentrations ranging from 10 mg/ml to 0.1 mg/ml. The concentrations are simulated as a weighted mixture of the lung attenuation material and pure molybdenum, to achieve the given concentrations. The range of concentrations is chosen to cover the experimentally observed lung accumulation of 1.5 mg/ml molybdenum NPs [1]. We note that XRF emission from the lungs experience severe self-absorption in the surrounding ribs, making self-absorption correction highly important. XRF-GPU is then used to simulate XRF tomography of the mouse phantom. Table 1 shows the acquisition parameters, typical for the laboratory arrangement described in Sect. 2.2.2. Note that short exposure times are used here, resulting in realistic scan times and very low radiation dose (c.f., Table 1). The dose is determined by storing the energy deposited by the absorbed photons in each voxel throughout the whole simulation.
The acquired simulated data is then reconstructed using our dual CT & XRF reconstruction algorithms, described in Sect. 2.2.3. The lung is segmented from the XRF reconstruction using the lung from the original mouse phantom as a mask to provide a histogram over the reconstructed concentration inside the whole 3D lung. Since the mean value mean (c ) of the distribution does not capture its shape, we introduce two additional metrics to quantify the accuracy of the reconstructions. First, we introduce the peak value of the histogram peak (c ) , and secondly its full width at half maximum FWHM (c ) . Together with mean c , these provide valuable insight into the shape of the histograms and hence, the quantitative uncertainty of the reconstructions. The mean and characteristic peak values can then be used to define error metrics relative to simulated concentrations sim (c ) , such as the

Quantitative tumor simulations
Localizing tumors using NPs is a major driving force in the development of small-animal XRF tomography [1,6]. Here, we investigate the performance of our laboratory arrangement to image few-mm sized spherical tumors surrounded by highly absorbing bones and assess the quantitative accuracy of the XRF reconstruction in the complex realistic environment. We modify the original mouse phantom and create a tumor-bearing model of the mouse by inserting three spherical tumors centered in the same axial slice in the hip region. Two sizes of tumors were investigated: a larger set of 3-mm diameter spheres, and a smaller set of 1-mm diameter spheres (both at 100 µm resolution). Figure 3 shows the tumor-bearing mouse phantoms. The simulated Mo NP concentrations in the tumor was chosen to 2 mg/ml, which is in the range of experimental observations, 2-7 mg/ml using passive strategies [9,31] and 3.5 mg/ml using targeted affinity ligands [32]. In addition, we simulate a very low concentration of 0.1 mg/ml for the larger tumors. The same acquisition settings (cf. Table 1) were used as in the lung simulations, with one difference: the slice separation was 200 µm instead of 500 µm to achieve better vertical resolution, owing to the small size of the tumors.
As the tumors are placed in different positions in the hip region, some will inevitably experience more severe self-absorption than the others. This poses the challenge in correctly quantifying all tumors, which are all artificially "injected" with the same concentration. Here we demonstrate the importance of correcting for self-absorption by reconstructing the data sets with and without correction.

Detecting sub-mm tumors
Early tumor detection is important for a number of reasons, one being that it increases the treatment prognosis. Ideally, pencil-beam laboratory XRF tomography should provide the resolution needed to observe sub-mm features such as early stage tumors. Sub-mm feature detection has been demonstrated experimentally with our laboratory arrangement before [1]. Here, we investigate the limits of detecting sub-mm tumors with our laboratory arrangement. We shrink the tumors used in the previous section to roughly 200 200 200 x x µm 3 volume each. The exposure times are increased to 100 ms to allow the small features to be detected, then we investigate at which concentration all sub-mm tumors can be detected. We note that the increased dose of ~300 mGy (c.f., Table 1) is still within the acceptable range of smallanimal CT [33]. For these small features, we value observability over correctly quantifying their concentration. We evaluate observability by two means; locally according to the Rose criterion [34] (stating that a squared signal-to-noise ratio ( 2 SNR ) of 25 > is required for a feature to be distinguishable by the human eye) and globally through the possibility of putting a lower threshold on the whole reconstructed volume to distinguish the true features from background noise.  Figure 4 demonstrates the effect of self-absorption and our correction method. The differences can be observed both qualitatively and quantitatively, here demonstrated on the 2 mg/ml Mo NPs simulated lung data set (c.f., acquisition parameters in Table 1). As seen in Fig. 4(a), not correcting for self-absorption leads to the incorrect conclusion that the core of the lung has accumulated much lower concentration of NPs than the outer parts. Furthermore, the outer parts are incorrectly quantified as well, as they still suffer from some self-absorption (although less than the central parts). This is quantitatively depicted in the corresponding histogram in Fig. 4(b), leading to lower reconstructed concentrations overall in the whole lung. Similarly, we show an example of potential over-estimation of the self-absorption, e.g., due to incorrect assumptions of various parameters, resulting in too much correction. The latter leads to the core of the lung falsely exhibiting higher accumulated concentrations than the outer parts. If the right correction is implemented, the homogenous distribution should be recovered, ideally leading to a narrow distribution around the true value (c.f., Fig. 4(b)). We note that these simulations can therefore be used to validate the physics model in the selfabsorption correction and to ensure that no over-estimation is made. We see that our correction attempt leads to a distribution with its peak near the true value, meaning that a large portion of the voxels inside the lung have a value close to the true value. Still, some voxels have significantly lower or higher concentrations than the true value. We note that in addition to the self-absorption correction itself, there are several factors that influence the reconstruction quality, such as exposure time per pixel and number of rotations in the tomographic acquisition. Nevertheless, the results indicate that there might be room for improvement in our correction method and that new attempts can be compared and evaluated easily using the present XRF-GPU computational model.  ± for 10 mg/ml, and 20% ± at 4 mg/ml. However, the mean values do not coincide with the peak values, and are typically lower for high concentrations. This means that the distributions are shifted towards lower concentrations for these data sets.

Quantitative large-organ simulations
For the lower simulated concentrations, the relative width of the distribution peak increases, and at the lowest concentrations the peak metrics seem to lose meaning. This can be observed from Fig. 5(b) which is the histogram for the lowest concentration (0.1 mg/ml). Clearly, the shape of the distribution has changed drastically compared to what is observed at higher concentrations (c.f., Fig. 4(b)). In other words, if one were to use the peak of the distribution as an indicator of the true concentration inside the lung, one would arrive at incorrect conclusions. While mean c is only 20% lower than the expected concentration, it is safe to say that the quantitative ability of our reconstruction algorithms breaks down at these low concentrations given the acquisition settings in Table 1.
The simulation results provide some insight into the behavior of our reconstruction algorithms. The results indicate that at high concentrations, peak c should be used as an indication of the accumulated concentration. At very low concentrations one should instead use mean c as an indicator of the true concentration, even though the reconstruction exhibits large variations. We note that while this shows that our reconstruction algorithms have room for improvement, this kind of insight would be impossible to achieve if not for the highly realistic simulations made available by the XRF-GPU computational model presented here. antification stud reconstruct fe y. Furthermore how importan rection and Fi y more bones, than the other g the histogram hey have limi orrectly within higher degree metrics presented hibiting a different dy. Figure 6 sh ew-mm sized t e, the shape an nt the correctio ig. 6(b) with c can falsely be r tumors if cor ms of each tum ited overlap. W n 5-10% relati (c.f., Fig. 6 hows that tumors at nd size of n of selfcorrection e believed rrection is mor in Fig. When the ive mean )). The bottom row shows the corresponding tumor histograms. Figure 7 shows the results from the same tumors but with significantly lower Mo concentration (0.1 mg/ml). Here, again we observe the importance of self-absorption correction, significantly reducing the RME. As observed earlier in the lung simulations, at low simulated concentrations the distribution tends to lower values, although this is somewhat counteracted by the self-absorption correction.
Distinguishing the tumors from the background noise is also facilitated through the selfabsorption correction (c.f., 3D-visualizations in Fig. 7). For self-absorption corrected data ( Fig. 7(b)), a higher threshold can be set on the XRF reconstruction to remove background noise (false-positives) without losing any significant tumor feature. Nevertheless, the tumor demarcation is naturally less distinct in comparison to the 2 mg/ml tumors (c.f., Fig. 6), both regarding size and shape. We note that this is due to a combination of low concentration and short exposure times, i.e., low dose.  Figure 8 shows the results from the 1-mm diameter tumors, artificially "injected" with the realistic concentration of 2 mg/ml Mo NP. As the tumor volumes are small compared to the 200 µm 3 voxel size of the reconstructions, the histograms provide limited information. Nevertheless, we again observe the trend that self-absorption correction significantly improves the average mean values in each tumor. It is interesting to note that the relative mean error (RME) in the estimates are within ± 10%, which agrees very well with the same concentration in the much larger lung volume (c.f., Fig. 5(a)). Furthermore, the tumor morphology in the reconstructions agree well with the simulated model.
The results in this study indicate that it is primarily the NP concentration per unit volume that determines the quantitative accuracy in our reconstructions. This is supported by the fact that the RME of the reconstructions corresponds well for both tumors and lungs at 2 mg/ml. Similarly, for the low concentration of 0.1 mg/ml, the RME of the 3-mm tumors were in agreement with those of the lungs.  Figure 9 shows the reconstruction of 3 200 200 200 μm x x tumors artificially "injected" with 2 mg/ml Mo NPs. Self-absorption correction has been implemented in the reconstruction. While the tumors themselves should ideally fit in a single reconstruction voxel (isotropic acquisition step size of 200 μm ), the XRF signal from the tumors is distributed among adjacent voxels (in a 3 by 3 by 3 voxel region) due to the nature of acquiring data with a step size equal to the feature size studied. In order to quantify the local 2 SNR , we integrate these regions to find the mean signal of each tumor and quantify the standard deviation of the background locally by studying a larger 9 by 9 by 9 voxel region around the center of each tumor (excluding the tumor itself). In addition to the 2 mg/ml tumors, tumors with 1 mg/ml Mo NPs were simulated. Table 2 shows the calculated local 2 SNR of the 1 mg/ml and 2 mg/ml tumors. Here we see that by the Rose criterion ( 2 SNR 25 > ) [34], the 2 mg/ml tumors should all be well observable. Due to its surrounding, the center tumor experiences a higher degree of selfabsorption, explaining the lower 2 SNR . In addition, Fig. 9 shows that these tumors can be completely distinguished globally from background noise by setting a lower threshold in the XRF reconstruction. In other words, it is possible to eliminate "false-positives" and still observe the true features. This is not true for the 1 mg/ml tumors. While the local 2 SNR of the center tumor is on the limit of observability (c.f., Table 2), it was not possible to set a background threshold without losing the feature itself or introducing false-positives from background noise. In other words, the center tumor artificially "injected" with 1 mg/ml could not be distinguished from the background. Overall, we note that these results agree with previous experimental results [1].

Sub-mm tumor detection
While these results should not be considered absolute, as they are highly dependent on the acquisition settings simulated, they demonstrate that there is a limit in observing small low contrast features in a complex and highly absorbing surrounding. Here we demonstrate that these limits can be realistically investigated with our XRF-GPU computational model. Fig. 9. 3D-rendered overlay of CT (grey) and XRF (red) reconstructions of sub-mm tumors artificially "injected" with 2 mg/ml Mo NPs. Self-absorption correction has been implemented. The inset zooms in on the left tumor, showing the sub-mm size in the reconstruction. Reconstructions are visualized with linear interpolation for facilitated qualitative inspection. The lower threshold in the XRF reconstruction is be set to completely distinguish the tumors from the background noise.

Conclusion and outlook
In this paper we have presented a realistic small-animal GPU-based computational model for XRF tomography. The core of the model is our very fast Monte Carlo simulation tool (XRF-GPU) together with its capability of simulating high-resolution realistic small-animal phantoms. The XRF-GPU computational model can be used to simulate any acquisition geometry and allows for a highly flexible geometry design including x-ray source and detector specifications. As it is highly accelerated, tomographic simulations are made possible in reasonable time-frames without high demands on computational resources. While a laboratory system for XRF tomography can be used to do experimental smallanimal imaging, evaluating the quantitative accuracy of the reconstructed data directly is practically impossible, since the true bio-distribution of NPs is unknown. This is where our presented model can be used to obtain a priori information of the quantitative uncertainty in reconstructions and predict observable features in a realistic setting.
To demonstrate the XRF-GPU computational model, we use it in a case study to simulate our own laboratory system. This way, we predict the quantitative performance of our reconstruction algorithms in a small-animal imaging scenario. This gives us key insights into its behavior in reconstructing different NP concentrations and shapes. We noted that this type of evaluation would not be possible using neither experimental phantoms (as they will always be oversimplified) nor with other simulation tools available today (too slow). Furthermore, we demonstrated the importance of correcting for self-absorption of XRF to achieve quantitatively correct results in realistic small-animal imaging settings. This is especially important for NPs based on elements with XRF energies experiencing significant attenuation in a few centimeters of soft-tissue, such as molybdenum.
We learn from the case study that accumulated concentration per unit voxel seem to be a more important factor than feature size in affecting the quantitative accuracy of our reconstructions. Furthermore, we demonstrate that that our system is capable of detecting submm tumors in the 200 µm-range artificially "injected" with realistic Mo NP concentrations, completely distinguished from background noise.
The results have shown that the physics model in our reconstruction algorithms can be further developed, to deal with the statistical nature of the low-concentration problem where only few XRF photons are detected in each projection pixel. Hopefully, this should lead to reduced uncertainty in reconstructing lower concentrations of NPs overall.
We note that while the XRF-GPU computational model is accurate in itself, the biodistributions simulated in the case studies are not fully realistic. In reality, tumor localization is more difficult since NPs typically accumulate in surrounding healthy tissue as well. Furthermore, while homogeneous distributions are practical for these types of evaluations, they are highly unlikely to be found in a real imaging scenario. However, more complex biodistributions can easily be artificially created and simulated to gain more insight into the performance of the imaging system. Finally, we note that our XRF-GPU computational model can be extended also to digital human phantoms. We predict that this tool can provide key insights into the design of optimal imaging geometries for XRF tomography of small-animals as well as humans.

Funding
The Wallenberg Foundation.