An image formation model for Secondary Ion Mass Spectrometry imaging of biological tissue samples

Secondary Ion Mass Spectrometry (SIMS) can provide distribution images of elements and molecular fragments with high sensitivity and spatial resolution. This study aims to exploit the potential of this modality as an imaging technique for biomedical applications. A model of image generation was developed and validated on experimental SIMS images. The model allowed for the selection of standard distance deviation (SDD) and nearest neighbor index (NNI) as suitable indices for the characterization of SIMS images, as they have been associated with sample morphology. Two regression models were proposed to correlate the SDD index and NNI with an index of effectiveness and acquisition parameters. The SDD index, due to its linear relationship with the image noise parameter, was less sensitive to noise. The model was then applied to study the effect of instrumental and analytical parameters, such as pre-sputtering time, on


Introduction
Secondary Ion Mass Spectrometry (SIMS) is a unique chemical imaging modality that is able to detect single isotopes and molecular fragments from a solid sample surface, whether inorganic, polymeric, or biological, combining chemical and spatial information to generate very high spatial resolution. SIMS is a hyper-spectral imaging technique since each image pixel consists of a complete mass spectrum; this feature allows the collection of a large amount of information by a single map acquisition but involves challenges in data handling and interpretation.
Although SIMS was introduced in 1962 [1], its first biological applications appeared several years later [2]. In 1989 the first study to provide a subcellular quantification method for SIMS ion microscopy images was published [3]. Later, the SIMS methodology became an active area of research, including cellular imaging [4], elemental and stable isotopic tracer distributions [5], characterization of boron distribution in Boron Neutron Capture Therapy [6][7][8], and investigation of calcium flow and stores in cancer research [9,10].
Analysis methods for SIMS data were proposed in 2006 to allow for mass spectrum processing by means of multivariate analysis methods [11][12][13][14][15] in order to correlate a large number of mass signals from each image. In these methods the spatial distribution of the SIMS signal was not taken into account as the whole data set was processed. Methods for quantitatively processing chemical images were proposed in [16], while an accurate mass scale calibration procedure was investigated in [17]. Recently, Time-of-flight Secondary Ion Mass Spectrometry (ToF-SIMS) [18] was employed to characterize cellular matrices and membrane lipid distribution [19][20][21][22], examine the immobilization process of PNA (Peptide Nucleic Acid) and its ability to hybridize to unlabeled complementary DNA fragments, and image elements and molecules in complex biological samples [23,24]. In particular, ToF-SIMS has been applied to protein evaluation on biodevices using the mutual information theory [25].
All of these previously described approaches were mainly focused on the spectral content of SIMS data rather than on the spatial information also provided by SIMS. In the present study, a model of SIMS image generation was developed. Using this model, we define indices capable of characterizing the elemental spatial distribution. We validated the model using SIMS images of biological samples.

SIMS technique and image generation
The diagram in Fig. 1 shows the typical components of SIMS instrumentation and the main analysis results. SIMS utilizes a partially destructive bombardment (sputtering) of the sample surface, in high vacuum conditions, by means of a focused primary ion beam (Fig. 1a) that sequentially scans the surface area of the sample (Fig. 1b) as a microscope "probe" in accordance with a discrete twodimensional pattern. The extraction of the species is performed by the subsequent production of secondary ion beams (Fig. 1c) according to a typical pattern representative of the chemical composition of the specimen surface. Secondary ions can be detected by two main approaches: Time-of-flight (ToF) or a magnetic sector mass analyzer. To image molecular fragments in biological samples, Time-of-flight Secondary Ion Mass Spectrometry (ToF-SIMS) is the technique of choice and will therefore be our focus in this study. In ToF-SIMS analysis, secondary ions are "gated" from the source region by an electrical field pulse and accelerated down the ToF flight tube (Fig. 1d). Low [m/z] ions travel at a higher velocity and reach the detector quicker than the slower ions with a high [m/z]. Calibration of the accelerating field and resulting flight times permits mass analysis for unknowns. In particular, the secondary ions are accelerated in the mass analyzer and the time t elapsed between the sample surface extraction and an ion reaching the mass analyzer is measured according to: where L represents a characteristic spectrometer geometric dimension, V represents the acceleration potential, and the product q = ez represents the charge of the ion. As the ions all possess the same kinetic energy, each ion acquires a characteristic velocity that depends on its mass m and thus reaches the detector at different times (t). In order to obtain a spectrum of the intensity vs the massto-charge ratio, all secondary ions extracted in a single excitation process are sequentially acquired for each pixel. The production of distribution images, based on mass spectra, makes SIMS an ion microscopy modality. With proper calibration, the SIMS imaging technique is capable of providing relative images as well as a semi-quantitative measure of the distribution of elements and molecular fragments with a spatial resolution of about 1 m, a sensitivity of up to g/g, and a mass resolution ranging from one hundred to one thousand atomic mass units (amu). The ability to provide chemical and morphological information at the same time makes SIMS a potential tool for analyzing the microlocalization of elements and molecular fragments in biological matrices, even at the sub-cellular level.
Due to the parallel detection of the mass analyzer, the complete mass spectrum is collected for each pixel (Fig. 1f). Thus, SIMS images contain both the mass composition and its distribution on the sample surface for each chemical species. Fig. 2 shows a typical example of a SIMS image from a biological sample.
ToF-SIMS can operate in both static (analysis is performed on the first molecular layer of the sample surface) and dynamic (performing depth profile scanning at a selected point of the sample) modalities. Static SIMS is limited to the analysis of the first monolayer. Because the primary ion dose is kept below 1E + 12 ions/cm 2 , the mass spectrum can reveal molecular information. Dynamic SIMS allows 3D analysis, exploring the entire sample depth. However, the primary ion dose may exceed 1E + 12 ions/cm 2 and only elemental and isotopic information can be obtained from the mass spectrum. The static modality (exploited in the present study) allows the acquisition of spectra by scanning a selected area of the sample with a selected array of discrete points (pixels) of approximately 2 nm on the first molecular layer of the sample surface.
Ion sources commonly used in ToF-SIMS are Ga+ or In+ liquid metal ion guns (LMIG). With these monoatomic primary species the imaging of biological samples was generally limited to low mass ions such as the phosphocholine head group at m/z 184, due to the large fragmentation of biological molecules, like proteins, lipids DNA/RNA. Cluster ion sources have the potential to increase secondary ion yields by up to 1000-fold in some cases.

Point pattern distribution analysis
A SIMS image is inferred from a counting operation, which is in contrast to most medical gray-level images that are obtained by means of a quantization of a continuous physical phenomenon. As shown in Fig. 2, a SIMS image could be modeled by a lowintensity two-dimensional point distribution. Hence, it could be  useful to introduce some analytical instruments borrowed from the so-called "Point Pattern Distribution Analysis" [26][27][28][29][30][31][32][33][34][35][36][37], which was designed to analyze this kind of point distribution.
One class of indices used in the literature for point distributions includes the so-called "Centrographic Statistics", which substantially account for the centrality and the dispersion of a spatial distribution. In particular, the standard distance deviation (SDD) represents the two-dimensional generalization of the standard deviation of a single variable in accordance with the following equation: where (x i , y i ) and (x c , y c ) represent the Cartesian coordinates of the i-th point and of the mean center value of the N point distribution, respectively. N represents the number of samples. The Shannon entropy index accounts for the gray level variation modalities in the histogram of an image.
As far as point pattern distribution analysis is concerned, two categories of methods can be recognized in the ecology and astronomy literature: the area-based (or Point Density-based) methods and the distance-based (or Point Interaction-based) methods. The former includes "Quadrat analysis" [29] and other techniques that consider the frequency distribution of the number of points belonging to a regular sub-region of the study area. The latter category contains the methods that characterize a spatial pattern by means of the relative distances between points and specifically includes nearest neighbor analysis [26][27][28], which involves the calculation of the nearest neighbor index (NNI). A generalization of the nearest neighbor analysis technique was proposed in 1976 by Ripley [33] and has been applied in many recent papers [30][31][32][34][35][36][37]. This is a multi-scale method that computes a function (or an index) called Ripley's K in order to identify and locate patterns, clusters, and other local distribution heterogeneities.

SIMS image analysis
A fundamental question in the analysis of SIMS images from biological samples if the detection of preferential patterns of deposition for a particular chemical specie, patterns that could be hidden by several noise sources. If the number of counts is high enough, the problem can be solved using the standard image processing techniques commonly adopted in medical imaging. This is the case of the SIMS image in Fig. 2b and c. Fig. 2c represents a hole in the sample (the silicon in the substrate is well visible). However, several species of interest may be characterized by a low number of counts. An example is provided in Fig. 2d that describes the distribution of the 99m Tc-NOET traces (a radioactive tracer commonly used in Positron Emission Tomography studies). A lower number of counts can be expected in the tissue hole (Fig. 2e), however the visual detection of this pattern is not trivial if the hole location is unknown (Fig. 2d). Point pattern analysis may solve this issue, detecting by an objective measure if a count distribution contains a significant pattern. It could be accomplished by comparing the observed pattern with random pattern generated by an appropriate model. Moreover, point pattern distribution analysis may be useful in the planning of a SIMS experiment, providing the number of counts needed for the detection of a pattern of a defined size, and understanding how image preparation method (as pre-sputtering) could affect the final results. In the following, a SIMS image model is developed and validated and several point pattern distribution indices are evaluated. Two examples of possible application of the developed method are proposed.

Image model
In order to describe the SIMS image generation process, a model was developed (Fig. 3) and implemented in the IDL 6.1 software environment. As previously described, the value of each pixel in a SIMS image of a metabolite represents the number of counts for a specific spatial location. The model assumes that the acquired SIMS image for a spectral peak, denoted as I, can be obtained by the application of a g h,1− function to an ideal pattern I 0 , as shown in Eq. (3): The ideal pattern I 0 consists of piecewise constant regions, each characterized by a value representing the concentration of a metabolite in the region. Hence, I 0 represents the "true" distribution of metabolites on the sample surface.
The g h,1− function accounts for the acquisition physics that depends on the instrumentation, and is independent of the selected spectral peak. In other words, g h,1− is the same for all the imaged chemical species. The g h,1− function includes a smoothing operation of the ideal pattern, which is dependent on a bandwidth h parameter, to simulate the uncontrolled ion redistribution and displacement due to the primary ion beam sputtering action as well as a normalization to the maximum value in order to obtain a probability map of count localization in the image surrogate space.
Kernel smoothing, a non-parametric regression method introduced by Nadaraya and Watson in 1964 [38,39], was used to implement the smoothing operation of the g h,1− function. Thus, the smoothed image I 0 is obtained as: where N represents the point number of the spatial distribution (i.e. non-zero image locations), the u vector contains the spatial coordinates of the point where the function is computed, X represents a N × 2 matrix including the spatial coordinates of all image locations, Y represents a 1 × N matrix including the intensity value in each location (i.e. the gray level in each pixel), K(v) is the kernel function, h is the smoothing parameter and a assumes a constant value of 2/ in a two-dimensional problem; in a two-dimensional problem, the X[i] notation indicates the two-component vector of the i-th point coordinates, extracted from the X matrix. The smoothing parameter (or bandwidth) h controls the smoothing degree as it changes the weight function (or kernel function K(v)) width. Gaussian and parabolic forms [40] have been evaluated for the kernel function. I 0 is normalized between 0 and 1 and used as a probability map to generate the g h,1− function result (Fig. 3). Hence, defining q as the total number of counts of the final surrogate image (1 − )q counts are generated with a probability to be placed in an image location proportional to the value of the normalized map in that location.
In order to identify and locate eventual structural discontinuities in the image, we deduced an analytical derivation formula in closed form. An example of the matrix resulting from the application of the kernel smoothing theory is shown in Fig. 4.
The second term of Eq. (3) (denoted by N ) accounts for all kinds of deviation from the ideal pattern I 0 due to the peculiarities of various chemical species. It includes species mobility, artifacts introduced during sample preparation, and environmental contamination. Hence, N depends on the chemical species imaged as well as the sample analyzed. The N term is obtained by a realization of a process that generates points with randomly distributed coordinates on the whole image domain.
The parameter (constrained to be between 0 and 1) represents the percentage of the total intensity q of the surrogate SIMS image generated by the N component. Hence, the total intensity generated by the g h,1− term in Eq. (3) will be (1 − )q. The final surrogate image with total intensity q is obtained as a sum of (1 − )q and q intensities.

Comparison of the image model and the real SIMS images
The developed model can be compared with real SIMS images of biological samples. The availability of microscope images perfectly aligned with SIMS images allows the generation of an I 0 image that perfectly matches the real phantom shape. Using Eq. (3), surrogate images can be generated depending on model parameters and compared with the real SIMS images.
The parameters h and in the image formation model can be set in order to maximize the similarity between the generated surrogate image and the original with the total intensity of each being equal. To compare the original and the substituted image, a threshold mask technique can be used. A square kernel filter scanned the silicon image (since it is the main constituent of the support wafer), shifting a few pixels at each step. The total intensity in the window is compared with an experimental threshold to generate a binary mask that can then be used as a mask for all images extracted from the same mass spectrum. A double size window scans the original and the substitute images in order to classify each square portion in three zones, depending on the total, partial or null intersection of the window with the masked ones. In each position of this window, the mean value of the intensity is computed. The mean value in the three zones is averaged and used as an index of image similarity. This process can be repeated for several chemical species.

Sample preparation and imaging with static ToF-SIMS ion microscopy
Images of biological samples were acquired during an imaging protocol developed to investigate the influence of microvascular alterations on the onset and development of myocardial dysfunction.
The distribution and accumulation of two flow tracers were studied using the ex vivo beating heart model of Langendorff according to different protocols of perfusion in order to reproduce low-flow ischemic conditions [41]. In the first experiment, a 250-L infusion of 99m Tc-NOET tracer solution (0.17 mM) was administered over a period of 2 min in an isolated mouse heart, with an expected uptake of 80% based on previous results [42]. In the second experiment, gadolinium-diethylenetriaminepentaacetic acid (DTPA) was infused. The whole heart was then frozen in isopentane cooled to −150 • C in liquid nitrogen [43]. Afterwards, the sample was embedded in OCT (Tissue Tek), mounted in a cryomicrotome holder, and sliced into 8-m thick coronal sections. These sections were placed in silicon holders (silicon wafer, single side polished, thickness of 1 mm, P type conductor, Aldrich) and moved on dry ice to a vacuum freeze dryer where the sample remained for several hours. Slices of different portions of the heart were analyzed to map the distribution of the tracer using SIMS, and the data obtained was correlated with complementary morphological information from optical images obtained with SEM (Scanning Electron Microscope) or optical microscopy. ToF-SIMS was performed with a Physical Electronics PHI TRIFT III, which is configured with a pulsed Ga + liquid ion gun and pulsed liquid Au + , Au 2 + , and Au 3 + ion guns. The Ga+ ion gun was used as reference ion source in the study. The TRIFT III model operated in static mode, which allows for the analysis of both elemental and molecular species at the sample surface. A 300 m × 300 m area was imaged. The mass range was set to 10,000 amu. The tracer 99m Tc-NOET is a well-known deposit tracer that carries a Technetium-99 nucleus ( 99 Tc is a long-lived radioactive fission product with a half-life of 2.13 × 10 5 years) that is used as the "probe" for SIMS analysis, meaning that the 99 Tc and molecular fragment masses of the molecule were mapped onto the sample surface. The 256 × 256-pixel SIMS images (.raw format) were extracted from the mass spectrum using the WinCadence (Physical Electronics) software environment.

Results
The developed SIMS image model could be useful in various applications. In this section we first report the results of the validation of the image model introduced in Section 3.1 against real SIMS images acquired as described in Section 3.3. Image comparison was performed following the approach in Section 3.2. In Section  Fig. 5 shows a typical example of the procedure described in Section 3.2. In particular, Fig. 5i shows the I 0 images extracted from the sample by the optical image in Fig. 5j. Fig. 5a-h shows SIMS images of several species, i.e., total ion (a), silicon (b), two phosphatidylcholine fragments (C 5 H 12 N + , C 5 H 14 NO + ) (c,d), carbon (e), sodium (f), potassium (g) and the TcN fragment of the 99 Tc-NOET tracer substance (h).

Model validation
In all studies, the optimal smoothing parameter h value was found independent of the chemical species selected and was equal to 3.0. This confirms the hypothesis that this parameter represents the ion dispersion due to the sputtering process. In order to maximize the similarity between images, the parameter had to be changed in accordance with the chemical species and the sample considered. Optimal values were <0.01 for b, f and g chemical species. This means that for images with a large number of counts the weight of the "noise" term in Eq. (3) is negligible. The other values were = 0.8 for C 5 H 12 N + (c), = 0.25 for carbon (e), and = 0.2 for 99 Tc-NOET (h). Comparison of values for C 5 H 12 N + (c) and 99 Tc-NOET (h) revealed that "extrinsic" tracers such as 99 Tc-NOET were less sensitive to noise than species in biological tissues, that had the same number of counts (36 vs 29, respectively). Visual comparison of Fig. 5c and h confirmed this fact.

Index quality assessment
Confidence interval boundaries for indices described in Section 2.2 were evaluated by means of a surrogate data model that generates a large number of random realizations. The model in Eq. (3) was used to generate grid-shaped intensity images of 200, 500, 1000, and 5000 counts. Grid bars with different spacings (10, 20 and 40 pixels) and thicknesses (2 and 10 pixels) were used as shown in Fig. 6. The parameter h was kept constant (h = 3.0) as previously described. Images with the parameter ranging from 0 to 1.0 in 0.1 increments were generated. In this way, the ability of different characterization indices to detect the grid was evaluated depending on the grid geometry, the amount of noise, and image intensity.
For each model configuration, 1000 realizations were generated and the distribution (i.e., mean and SD) of the characterization indices was evaluated on surrogate images. The statistical distribution of indices was compared with the distribution computed from a "noise only" image (i.e., an image with a random point distribution with the same intensity as the surrogate image). Fig. 7 plots the highest noise levels that allow the proposed indices to detect grid geometry, as we expected that the noise level  would modulate the significance of an "effective" index. The most effective indices for the characterization of different grid geometries, and presumably of imaged object dimensions, were SDD and NNI. Thus, we subsequently focused on these two indices. Figs. 8 and 9 show the behavior of the SDD index and NNI behavior as related to model configuration. The two indices both lose statistical significance compared to the random distribution for a value that is a decreasing function of the grid bar thickness and an increasing function of grid bar spacing.
Two regression models to fit the SDD and NNI curves were identified in order to have a closed form analytical expression. Fig. 9. Nearest neighbor index (NNI) trends as related to noise level ( factor), grid bar thickness (first digit in label) and grid bar spacing (second digit in label). The two parallel lines indicate the upper and lower bounds of the 95% confidence interval of NNI for a "noise only" image. For the SDD curve, a linear regression model was determined, but for the NNI an exponential regression model achieved the best fit. Figs. 10 and 11 show the regression models inferred from several model geometries.

A model application: pre-sputtering simulation
In this section, we show how our model can be exploited to investigate the influence of pre-sputtering on image quality. Pre-sputtering is an important preprocessing step routinely performed in SIMS. It consists of a uniform bombardment of the first  molecular surface layer of the sample to remove any environmental contaminants that settled on the surface during its preparation. This operation is known to potentially degrade image quality by introducing noise.
Spectra of the same sample were acquired for different presputtering times (0, 60, 210 and 420 s, respectively). A current intensity of 2 × 10 −9 A in DC mode was used with a beam energy of 22 keV. The corresponding ion doses were 6.0 × 10 13 (60 s), 2.1 × 10 14 (210 s), and 4.2 × 10 14 (420 s). In particular, the ion maps of three phosphatidyl-choline fragments (shown together with the total ion map in Fig. 12) were analyzed. The Signal-to-Noise Ratio (SNR) for the images was computed as the ratio between the  total intensity in the tissue region and background region where assignment to the masked or unmasked regions was established by considering the silicon image as previously described. The total intensity trend, the decreasing trend of the SNR, and the increasing trend of the value with increasing time of pre-sputtering for both the total ion and the three fragment images are shown in Fig. 13. These trends agree with the hypothesis that the parameter represents image noise.

A model application: pattern detection
Finally, the previously introduced SDD index and the developed model are applied to analyze and detect a pattern structure within a mouse heart tissue. The aim is to verify the inhomogeneous distribution of the gadolinium-DTPA flow tracer, that is expected to diffuse only in the inter-cellular space. SIMS maps related to 50 main metabolites were extracted and the statistical significance of the SDD index was evaluated as previously described for each map. Fig. 14a and f shows the SIMS maps of gadolinium-DTPA and calcium, respectively. Only gadolinium-DTPA map revealed a statistically significant difference of SDD index (SDD = 37.57) in respect to the random distribution (confidence interval 32.04-36.33), demonstrating the ability of the proposed approach to identify "interesting" samples. SDD index for the calcium map was 35.46, with a confidence interval for the equivalent random distribution of 34.50-35.82.
The identified pattern could be better characterized by an image filtering procedure. A square filter of increasing kernel size was applied on the gadolinium-DTPA and calcium distribution maps and the SDD index was locally computed for each kernel size. Fig. 14b-e shows the SDD masks for gadolinium-DTPA image with kernel sizes of 5, 11, 15, and 19, respectively. The distribution pattern arises in Fig. 14d, confirming that gadolinium-DTPA diffuse only in inter-cellular space. Fig. 14g-l shows the SDD masks for calcium with the same kernel sizes. The uniform distribution of calcium is evident.

Discussion
This study attempts to introduce a new approach for the characterization of complex matrices, like biological tissues, with SIMS. A model able to describe the image generation process was developed and applied in few reference sample: the impact of the instrumental acquisition parameters, matrix effects and ions mobility on the pattern of distribution have been investigated and assessed using the index in the model. As shown during its validation (Section 4.1 and Fig. 5), the model is able to estimate the noise effect of different species during the acquisition process. This approach could be useful to forecast the pattern and extension of biological structures for a given chemical species according to a specific experimental set-up. Following this approach, taken the index of a chemical species/fragment selecting the relative mass number, the number of counts needed to obtain a significant difference in SIMS signals between two tissues of a given size can be derived.
Moreover as we presented in the paper, the standard distance deviation (SDD) index and nearest neighbor index (NNI) were shown to be able to describe the sample morphology. Hence, two regression models were proposed for the SDD index and NNI. The SDD index, due to its linear relationship with image noise, was shown to be more robust in a wide range of image quality conditions.
As a test of the model, we investigated the effect of presputtering time on SIMS image quality (Section 4.3). As shown in Fig. 13c, the model allowed the linear relationship between the presputtering time and the noise induced on the SIMS images to be shown. This approach could be useful in experiment planning, giving the optimal sputtering time to obtain readable images. Finally, we demonstrated the ability of the proposed approach in detecting significant patterns in a large number of metabolite maps and in characterizing the detected patterns.
This approach should be considered as complementary to other methods: variants of factor analysis (e.g. PCA, MCR and MAF) (Tyler, 2006, Tyler et al., 2007 have been developed to reduce the number of variables (i.e. ion peak areas in the case of ToF-SIMS) and are useful to handle complex data set and to identify relationships between variables and samples (i.e. image pixels). Some of these methods (i.e. PCA, MCR) could be sensitive to preprocessing and normalization and do not take into account spatial information. MAF is independent of preprocessing and normalization and includes information on the nearest neighbor of each pixel. The method here developed differs from these in several aspects. It involves the development of a simple mathematical model that mimics the physical phenomenon that origin SIMS images. As previously pointed out, the model allows to forecast the kind and size of biological structures that can be characterized. We also propose a new index for SIMS images (i.e. SDD index) able to describe the ion dispersion directly from the image counts. Hence, a quantitative comparison between our methodology and factor analysis is difficult due to the different information obtainable with the two approaches.
Some limitations of the method come from the sample preparation steps. Briefly our sample was obtained with 8-m-thick slices in order to increase the data quality and the spatial localization of the tracer: the thickness of the sample slices may produce unwanted charging effects and lead to a lower image quality. Cryomicrotoming may also produce significant alterations such as the "smearing" of the matrix sample during the cutting, and this also leads to a degradation of image quality; additionally the surface morphology of the sample can affect the generation of secondary ions. However, as previously mentioned, the two components of the model introduced in Eq. (3) cover several image degradation sources, such as limits in image resolution, matrix effects, sample contamination, preparation artifacts, etc. It may be of interest to develop a more detailed model that could take into account image degradation sources separately. Conversely, the simplicity of the model allowed us to effectively fit the model to real data for SIMS images related to chemical species with a low number of counts.

Conclusion
In conclusion, a model for SIMS image characterization complementary to well established approaches like MAF and PCAs was developed and applied to different biological samples. The model could exploit the potential of SIMS imaging for the investigation of the pattern and distribution of molecular fragments and bioactive molecules in biological tissues, allowing SIMS to be used in the pharmaceutical and biomedical fields, as well as in transitional research.