Improved single particle ICP-MS assessment using a novel Python-based data processing algorithm (Sparta) for nanoparticle quantification

Abstract

Single particle inductively coupled plasma-mass spectrometry (spICP-MS) is a valuable tool to characterise nanoparticles (NPs) regarding their element-specific mass, size and particle number concentration (PNC). However, spICP-MS still suffers from a lack of harmonised and transparent data processing algorithms, resulting in little user-flexibility in adapting parameters, when working with e.g. the manufacturer software. In this study, we present a transparent Python-based algorithm (called ‘Sparta’), validated and critically compared with existing data processing methods (SPCal and an in-house Excel method as well as two commercial instrument software), applied for measurements of ∼30 nm Au, ∼74 nm TiO₂ and ∼50, ∼100 and ∼300 nm SiO₂ NPs, using instruments from two different manufacturers using milli vs. microsecond dwell times. Sparta is capable of correcting baseline drift, determining the particle detection threshold (PDT) via the Poisson and iterative Gaussian method, performing a peak summation necessary for microsecond dwell times, and even extracting specific mass or size distributions from e.g. polydisperse materials via a Gaussian peak-fitting. Although all data processing methods benchmarked sizes and PNCs suit well for Au NPs, results show that millisecond dwell times systematically overestimated sizes for TiO₂ and SiO₂ (from 50–100 nm). For microsecond dwell times, only SiO₂ (50 nm) showed slight overestimation due to the methodological LOD_size of 53.1 nm for our algorithm. Nevertheless, Sparta accurately removes spurious background events of challenging samples such as SiO₂ at larger particle sizes (i.e., 300 nm). Thus, it can be readily applied to other engineered and natural NPs or even for biological cells (single cell ICP-MS) showing its great potential in improving data processing for spICP-MS.

---

Introduction

Single particle inductively coupled plasma-mass spectrometry (spICP-MS) is a widely used method to characterise engineered, inorganic nanoparticles (NPs) in aqueous suspensions. spICP-MS can quantify NP size and particle number concentration (PNC), having the unique advantage of breaking down single NPs to their elemental composition. To our knowledge, spICP-MS was first reported by Degueldre and Favarger, who presented the feasibility of analysing several colloids ranging from 150 to 400 nm.¹ Particularly, gold (Au) and silver (Ag) NPs reveal low size detection limits (LOD_size) of 12.3 and 9.0 nm,² respectively, using spICP-MS. Such low LODs are obtained due to the (i) known composition, size and morphology of these NPs, (ii) the few, well-known ICP-MS-based interferences affecting these elements and finally (iii) their low (ionic) background (i.e., low dissolution).^3–6 Also, more complex NPs, such as titanium dioxide (TiO₂) or silica (SiO₂), have been studied using this technique.^5,7,8 On the one hand, TiO₂ NPs have gained attention due to their use in food and especially in cosmetics such as sunscreen, consequently being released into the (aquatic) environment, posing a potential danger to living organisms.⁹ Investigating TiO₂ NPs is challenging as they tend to agglomerate in neutral, slightly acidic aqueous suspensions (zeta potential at pH 6.2: ∼0 mV), resulting in larger size distributions.¹⁰ However, suspensions may be bath sonicated (temperature-controlled to prevent dissolution) prior to analysis, to counteract agglomeration of NPs.^11–13 Furthermore, ICP-MS based interferences (e.g. ⁴⁸Ca⁺) on the most abundant isotope (⁴⁸Ti⁺) and lower sensitivity result in a higher LOD_size (∼30 nm) for TiO₂ compared to Au or Ag.³ On the other hand, SiO₂ NPs are also common additives for food and cosmetics to prevent products from clumping and binding together.^14–16 The analysis of SiO₂ NPs using spICP-MS is even more challenging than TiO₂ NPs due to the intrinsic high ionic and particulate background. This background originates from glass-containing equipment used in ICP-MS and the high number of polyatomic spectral interferences (e.g. ¹²C¹⁶O⁺ or ¹⁴N¹⁴N⁺) for the most abundant isotope (²⁸Si⁺).¹⁷ Overall, these challenges generally result in a poor LOD_size of 105 nm.³ It is important to enable the detection of even smaller NPs as the potential danger of all types of NPs for living organisms is not yet fully understood and it depends on their elemental composition and size.⁹ As cell incorporation is greater with smaller NPs, small NPs pose an especially high potential cytotoxic danger for living organisms, depending on their elemental composition. Therefore, there is a niche in research for more efficient methods to characterise and distinguish (small) NPs by their elemental composition that can be improved via spICP-MS.¹⁸

Although spICP-MS has several strengths such as high sensitivity in characterising single NPs and a large variety of information provided, some intrinsic challenges and limitations remain. For instance, to calculate NP sizes in spICP-MS, assumptions must be made regarding their composition (e.g. Ti as TiO₂ equivalents), density and shapes; NPs are typically assumed to be spherical. Moreover, quadrupole (Q)-based spICP-MS (used in this study) cannot predict the elemental composition in each single NP in contrast to time of flight (TOF) spICP-MS. spICP-Q-MS works sequentially, which means that the NPs in the sample are assumed to be homogeneous in composition and therefore contain either one element or a fixed ratio of different elements. Nevertheless, the latest quadrupole-based ICP-MS systems are generally ca. one order of magnitude more sensitive than TOF-based ICP-MS systems making spICP-Q-MS highly valuable for detecting small NPs (lower LOD_size).^19–21

Another challenge in spICP-MS is the need for a good interplay between dwell time (also known as integration time), i.e., the time during which data are collected, and sample dilution. When using dwell times in the millisecond range (most commonly 2–10 ms), even with a dilution to yield about 10⁴ and 10⁵ NP per mL,^11,22,23 more than one event can occur within a dwell time, leading to an overestimation of the NP size and underestimation of the PNC. A big step-forward during measurements was the introduction of microsecond dwell times, available with the new generations of ICP-MS.²⁴ As a typical NP duration, in the plasma until ionisation, is between 200 and 1000 µs,^3,25,26 depending on the NP size and elemental composition, dwell times between 50 and 100 µs (ref. 23) are typically preferred to separately detect NP events and allow higher PNCs, while maintaining the low occurrence of two-particle events and even decreasing the NP LOD_size.^3,24,27,28 A dwell time between, or slightly higher than, the NP duration should be avoided to prevent partial/split event detections.²⁴ As presented in another study for 60 nm Ag NPs, when ionic Ag (≥1 µg L⁻¹) is present, the distinction of NPs from the ionic background is only successful for dwell times ≤100 µs, as for millisecond dwell times, the background and NP histograms merge and prevent clear NP identifications.²⁴ Thus, using microsecond dwell time not only significantly reduces the impact of ionic content but also improves the precision from ∼5% (10 ms) to ∼1% (≤100 µs) and significantly lowers the occurrence of multiple NP events in an integration time.²⁴

However, spICP-MS analysis currently lacks of harmonised, standardised and, particularly, transparent data processing procedures.²⁷ The manufacturers of ICP-MS provide user-friendly software for simple case studies. However, the exact NP identification algorithm is often not as transparent (e.g., not clearly documented) and not all relevant settings can be freely adapted for more complex samples. The instrument software typically shows excellent results for well-defined, optimally mono-element NPs, as presented in an application note for 10 nm Au NPs (LOD_size 6.5 nm) with low (ionic) background/noise content.²⁹ However, when the (ionic) background increases, the frequency vs. signal histograms of the background and NPs merge and overlap. When analysing elements with a higher (ionic) background such as Ti or Si, the NP separation from the background is only successful when the ionic content and the ICP-based interferences are (very) low, as presented in application notes for TiO₂ (LOD_size ∼30 nm)³⁰ and SiO₂ (LOD_size ≤ 100 nm).³¹ As this is rarely the case, many laboratories develop in-house data-acquisition procedures to process spICP-MS data. Many users rely on their internally developed scripts to have the possibility of adapting parameters for the interpretation, based on their specific samples and conditions.²³ However, their detailed algorithms are rarely presented and seldomly available. There is a limited selection of free software available, such as SPCal,³² Nanocount^33–35 and calculation tools (e.g. RIKILT Wageningen UR), each has its own strengths and weaknesses but in specific systems limited in their adjustment options. spICP-MS data treatment approaches were critically evaluated in a recent study, where the authors summarised that spreadsheet-based software (e.g. Excel or Origin) have the advantage of a smooth learning curve, widespread use and provides quick results, but are limited in their capabilities regarding multi-element analysis and the difficulty of implementing iterative processes.³⁶ They further stated that open-source software (e.g. SPCal or Nanocount) are transparent and have more advanced treatment options than spreadsheets.³⁶ Though they are user-friendly for quick calculations, but in-depth understanding requires a steep learning curve and are often discontinued or not maintained.³⁶ These disadvantages are even more pronounced for programming options, as they usually have no user interface, a particularly steep-learning curve and are usually developed for special scenarios, but have the impregnable advantage of enabling complete control over all (advanced) treatment processes and theoretically no limitations remain.³⁶ In any case, currently, common data processing approaches for quantifying particle detection threshold (PDT) include averaging the whole dataset and collecting only data which are higher than the average (µ) + 3 or 5 standard deviations (SD) in an iterative way.^6,37 Other approaches such as SPCal implemented local data thresholding (so-called rolling median filter) to correct for baseline drift scenarios.³² This is frequently accompanied by deconvolution methods, as presented in a previous publication for 10, 15 and 30 nm Au NPs, to improve the background-free NP identification.³³ These deconvolution methods are usually based on Poisson and/or Gaussian fits.^17,32,33,37

In this study, we aim to provide our novel data processing ideas and improvements within our Python algorithm named ‘Single Particle Analysis & Reliable Tracking Algorithm – Sparta’ to the Scientific Community, freely available via GitHub (link in Table 1), explained and validated in this study for spICP-MS of engineered NPs. We present a transparent Python-based data processing algorithm developed for spICP-MS, equipped with quadrupole mass analyser(s), combining already established strengths and tackling remaining issues in a constructive way. Our algorithm is capable of distinguishing NPs from the background signal based on the conventional Gaussian method in an iterative way⁶ and via the Poisson method,^32,38 which is especially valuable for low-background elements such as Au. Furthermore, we present a baseline drift correction and a peak summation for both methods when using dwell times smaller than the NP duration combined with an outlier removal for NP artefacts such as agglomerates. Here, we compare and validate our novel algorithm with existing data-acquisition methods from two manufacturers (MassHunter from Agilent Technologies and Qtegra from Thermo Fisher Scientific), one freely-available software (SPCal³²), and an in-house algorithm based on Excel (Uni-Oviedo).^39,40 To increase the applicability of the algorithm and verify the instrumental performance, we further evaluate micro vs. millisecond dwell time analyses (100 µs vs. 5 ms) for the test samples (NPs of Au, TiO₂ and SiO₂ of different sizes).

Materials and methods

Materials and reagents

All ionic standards and particle suspensions were prepared in ≥18.2 MΩ cm ultrapure water (Barnstead GenPure UV-TOC, Thermo Fisher Scientific, Bremen, Germany). Ionic calibrations (0–50 µg L⁻¹) were prepared in 2 mL polypropylene (PP) microreaction tubes (Eppendorf, Hamburg, Germany) using ICP standards (1000 mg L⁻¹): Au as HAuCl₄ in 7% HNO₃, Ti as (NH₄)₂TiF₆ and Si as (NH₄)₂SiF₆ both in H₂O, pH ∼3.5, all from Certipur (Merck, Darmstadt, Germany).

Nanoparticle samples and preparation

The investigated particles varied between 30 and 300 nm in size. In this manuscript, we use the term ‘nanoparticles’ (NPs) for particles in the metric nano-range between 1 and 1000 nm. Three different types of certified reference materials (CRMs) were investigated: (1) gold (Au) NPs (citrate stabilised) with a nominal diameter of 30 nm and a modal diameter of 32.7 ± 2.0 nm (LGCQC5050, LGC Standards, Teddington, UK), certified for PNC and modal diameter via spICP-MS and particle tracking analysis, respectively; (2) titanium dioxide (TiO₂) nanopowder consisting of 76% anatase and 24% rutile (NIST SRM 1898, National Institute of Standards & Technology, Gaithersburg, USA), re-suspended in ultrapure water resulting in certified diameters of 71 ± 4 and 77 ± 7 nm via laser diffraction spectroscopy and X-ray disc centrifugation provided by the manufacturer, respectively, and PNCs calculated by the authors based on the certificate information; (3) aminated silica (SiO₂) nanospheres of 50, 100 and 300 nm, certified via transmission electron microscopy (TEM) and, suspended in ethanol (nanoComposix, San Diego, USA) for which both calculated PNCs and certified sizes are provided by the manufacturer's certificate. NP suspensions were diluted in ultrapure water to final concentrations of ∼5 × 10⁴ NPs per mL. For deagglomeration, all suspensions were sonicated (SONOCOOL 255.2, Bandelin, Berlin, Germany) for ≥15 min at 20 °C prior to spICP-MS measurements. An overview of the NP suspensions used, and their certified parameters are summarised in Table S1.

spICP-MS instrumentation

Two instruments were used for the benchmarking: the 8900 ICP-MS/MS (Agilent Technologies, Waldbronn, Germany) and Thermo iCAP-TQ (Thermo Fisher Scientific, Bremen, Germany), both equipped with a total consumption spray chamber and a microconcentric nebulizer (High Sensitivity Single-Cell Sample Introduction System for ICP-MS, Glass Expansion, Port Melbourne, Australia). For the sample introduction, a Fusion 100-X (Chemyx, Stafford, USA) or E-1000 (New Era Pump Systems Inc, Farmingdale, USA) syringe pump (Thermo or Agilent, respectively) equipped with a 1 mL syringe (Hamilton, Bonaduz, Switzerland) was used set to a constant flow rate of 10 µL min⁻¹. For the SiO₂ measurements with the Thermo instrument, a cyclonic spray chamber was used to reduce the glass-based background, equipped with a peristaltic pump with a constant flow of 438.1 µL min⁻¹. For comparison, these measurements were also repeated with the previously mentioned total consumption system, but the background was higher in this case, probably due to the additional use of the glass Hamilton syringe. Detailed information on the ICP-MS settings and tune parameters can be found in Table S2. In addition, the ionic limits of detection (LOD) and quantification (LOQ) using 3 and, respectively, 10 standard deviations (SD) + average (blank) criteria are presented in Table S3. Dwell times of 0.1 and 5 ms with total acquisition times of 40 s and 120 s, respectively, were used for all measurements using the Agilent instrument, whereas 5 ms dwell and 120 s total acquisition time were used for the Thermo instrument. Author's note: for the Thermo instrument (beta version) used in this work, the setting of short dwell times is a limitation, but short dwell times can also be selected for the new Thermo iCAP-TQ ICP-MS instruments. The transport efficiency (η), also known as nebulisation efficiency, providing the ratio of NPs arrived in the plasma divided by the NPs in suspension, was measured daily and individually determined for each dwell time and instrument using the 30 nm Au NPs via particle number method, as described elsewhere.^6,41

Comparison of methods and thresholds

The particle detection threshold (PDT) criterion and dwell times used in this study are shown for all validation and comparison data processing methods applied in Table 1.

Table 1 Data processing methods and their corresponding versions used for the interpretation of the spICP-MS data. The ‘criterion’ describes how the particle detection threshold (PDT) was set to identify particle events

Method/software	Version	Threshold method/criterion	Integration/dwell times (ms)
a ‘k’ is a factor commonly set between 3 and 7 (Sparta) but can be extended to any number if necessary. For SPCal, it can be adapted between 1 and 8.b For the Thermo instrument (beta version) used in this work, the setting of short dwell times is a limitation, but short dwell times can also be selected for the new Thermo iCAP TQ ICP-MS instruments.c The code of our Python algorithm (Sparta) (version 1.1) is available via GitHub: https://github.com/Steffen9204/Sparta-sp_scICP-MS.
Agilent, MassHunter 5.1 (measurement)/5.2 (data processing)	G7201D/D.01.02/Build 708.1	‘Automatic’	0.1, 5.0
Thermo Qtegra 2.10	2.10.3324.131	‘Automatic’	5.0^b
SPCal^32	1.4.5	Gaussian (µ + k^a SD) & Poisson; ‘automatic’ choice	0.1, 5.0
In-house algorithm based on Excel (Uni-Oviedo)^39,40	—	Gaussian (µ + k^a SD)	5.0
Sparta^c	1.1; Jupiter Lab (2.2.6) via Anaconda Navigator (2.6.0)	Gaussian (µ + k^a SD) & Poisson	0.1, 5.0

---

Data treatment for spICP-MS: in-house Python algorithm – single particle analysis & reliable tracking algorithm (Sparta)

In this subsection, the theoretical calculations used for Sparta (Gaussian and Poisson methods) are described and explained. Further calculations, which are not primarily used in this manuscript, are listed in eqn (S1) and (S2). All densities, mass fractions and ionic calibration parameters are listed in Table S4.

Baseline drift correction. In some cases, especially when the (ionic) background is high, a baseline drift is occurring (i.e., not constant over time) and the drift must be corrected to accurately account for the particle events. Based on the premise introduced in a previous publication,³⁴ we first apply a baseline correction to our raw data as shown in eqn (1),where I_corr. is the corrected signal, I_raw is the raw signal, Mode_{Iraw (±50)} refers to the moving mode of 100 of the closest signals and is the average of the detected moving modes of 100 of the closest signals. This is further expressed in Fig. 1 before further processing and calculating the PDT. Fig. 1a and b show the baseline drift of the raw data before the correction, which would complicate extracting particle events from the background (several local PDTs would be necessary). Fig. 1c and d demonstrate the corrected raw data, where particle events can now easily be extracted and separated from the background with only one PDT. The ‘moving mode’ represents the mode of each data point based on the 100 closest signals moving for each data point.

Fig. 1 Baseline-drift correction via Sparta for ¹⁹⁷Au⁺ (10 ppb, ionic standard in ultrapure water, 5 ms dwell time). (a) Raw signals (signal vs. time) before correction. (b) Raw histograms (frequency vs. signal) before correction. (c) Corrected signals (signal vs. time). (d) Corrected histograms (frequency vs. signal).

Particle detection threshold (PDT). Independent of the chosen dwell time, first, a PDT must be set. Two methods were implemented: the (i) Gaussian and (ii) the Poisson method. The detailed iterative Gaussian method (average (µ) + factor (k) × SD) is described elsewhere.⁶ Briefly, the SD and average were calculated from the entire dataset, and all data points above µ + k × SD were iteratively removed until there were none left above this threshold, and detected as particle events. Data points ≤ µ + k × SD were detected as background. Based on the recommendation from Laborda et al. (2019) to also consider 5 × SD instead of 3 × SD, to reduce false event detections, we expanded our approach for k values between 3 and 7.³⁷ For a not purely Gaussian background emerging from a mixture of dissolved elements but also particles (e.g. SiO₂ from glass), a higher factor k must be selected to maintain the same false positive rate, although a higher k value increases the risk of small particle events being overlooked. We used a higher ‘buffer value’ (max. 40 iterations) to be sure the iterations were sufficient as more iterations do not change the final threshold and the code operation is not time-consuming, even though we did not detect more iterations than 17, thus number ‘40’ was set as default to finally obtain the PDT for identifying particle events. In version 1.1 of Sparta, there is a possibility to ‘enable’ the by default disabled ‘clear_cut’ function if there is a clear separation in the histogram between the background and the NP signal distribution (for well defined, monodisperse distributions). If enabled, the iteration is stopped and the ‘cps value’, at which the frequency has reached ‘0’, was selected as the PDT to avoid disturbing background signals in the final particle output, as is often the case for low background elements such as Au. For the datasets presented in this manuscript, this function was, however, not used. For the Poisson method, we used the statistical basis from Currie (1968) for paired observations (alpha = 0.05 was assumed) and Lockwood et al. (2021), but expanded it by adding the average from all raw data before separation, as we often observed an underestimation of the L_C and PDT_Poissonfor our CRMs to have base accounting better for varying background, as presented in eqn (2) (adapted from ref. 32, 38, 42 and 43),and eqn (3) (adapted from ref. 32, 38, 42 and 43),where µ_raw refers to the average of all raw signals, ε is the correction factor, which is ‘0.5’ in case the average of all raw signals is <5 counts, otherwise it is ‘0’ and τ is the chosen dwell time (their use is explained in the ‘Peak summation’ section).^32,38 Both Gaussian and Poisson distributions are probability distributions, each fitting in specific cases better to the samples measured and thus both having their own value. As explained in a previous publication, quadrupole-based ICP-MS instruments equipped with electron multipliers follow a Poisson distributed noise for low-count signals, while Gaussian distribution often better suits high-count signals.³⁷

Peak summation. When choosing dwell times below the NP duration (i.e., 0.1 ms), a peak summation is necessary for both methods (Gaussian and Poisson) (Fig. 2a and b). Our algorithm performs a peak summation if dwell times are <2 ms (based on the typical NP duration). For the Gaussian method, a particle event is detected when a signal exceeds the PDT (green vertical lines) and is merged until a data point falls below the PDT. All signals below the PDT are recognised as background (grey vertical lines). The average of the background data is subtracted from each particle event (from each green vertical line, not from the final summed-up data) to set the baseline to ‘0’. This step is the reason the baseline drift correction is crucial. For the Poisson method, the event detection works similarly. The only difference is in the use of the critical value (eqn (2)). When a particle event exceeds the PDT and is thus detected, all data points that precede and/or follow it and that are already above the critical value are summed. That is, all the data points before and after the detected signal(s) that are above the critical value, will be summed until the data points are below the critical value and thus recognised as background (left peak in Fig. 2b). Events between the critical value and PDT are not detected as particle events and are merged to the background data (right peak in Fig. 2b). The advantage of choosing short dwell times (i.e., 0.1 ms) is that higher PNCs can be measured as shown in Fig. 2a. For the longer dwell times (i.e., 5 ms), only one event must occur in a 5 ms time to avoid double events, which would result in multiple signal intensities and thus spurious ‘larger’ NPs.⁴⁴ Thus, for the cases of 5 ms, events are accounted individually but not summed-up when appearing adjacent over time.

Fig. 2 Peak summation principle of Sparta for dwell times ≤ NP duration (i.e., 0.1 ms). (a) Gaussian (b) Poisson.

Outlier removal. For well-known systems, such as our model NPs (CRMs), having monodisperse distributions (as presented and used in this manuscript), the ‘Outlier removal’ was enabled to remove large agglomerates or other artefacts. However, we suggest for unknown, natural samples to disable (set by default) the ‘Outlier removal’ (available from Version 1.1 of Sparta), as detected ‘agglomerates’ could indeed correspond to large particles. If enabled, particle data above PDT are detected as outliers if they are > µ + i SD. This outlier filtering step is designed to remove only certain events that may distort some measures of central tendency, such as the average. The outlier detection is monitored using histograms. If enabled, the factor ‘i’ is set as default to ‘3’ but may be increased for well-known samples to any number if a high background is expected, such as for Si. One example where ‘i’ had to be increased to ‘12’ is shown for 300 nm SiO₂, 0.1 ms dwell time in Fig. S1.

Transport efficiency. The transport efficiency (η) was determined for both dwell times (0.1 and 5 ms) individually using three replicates each of the 30 nm Au NPs via the particle number method as shown in eqn (4) (adapted from ref. 6, 32 and 41),where DF refers to the dilution factor; T is the total acquisition time; V is the sample volume flow and PNC refers to the particle number concentration. Based on the background, the ‘threshold’ was daily adapted usually varying between 1000 and 5000 counts per second (cps). For instance, for a 5 ms dwell time (Agilent), a threshold of 3000 cps was chosen. The final η is the average of three replicates, which are listed in Table S5.

Particle number concentration. The particle number concentration (PNC) was calculated after the peak summation for both dwell times, either considering all events above PDT or only extracting chosen events under the Gaussian peak-fitting as presented in eqn (5) (adapted from ref. 32),

Particle masses and sizes. Masses emerge from the particle raw signals and represent the elemental mass per particle; sizes are estimated using assumed mass fractions and densities presented as size equivalents (e.g. Si as SiO₂). The calculation of the elemental particle masses is demonstrated in eqn (6) (adapted from ref. 32),and the calculation of the spherical particle sizes as expected size equivalents in eqn (7) (adapted from ref. 32),where m_elemental refers to the elemental particle mass; f is the mass fraction (e.g. for Si as SiO₂ equivalent: ); M is the molar mass and ρ describes the particle density.

Limit of detection (LOD) for masses and sizes. The LODs were calculated from the background signals, singled out during particle event selection (i.e., grey signals, Fig. 2a and b), in a similar way as for particle masses and sizes, just using instead of the particle raw signals the corrected PDTs (Gaussian or Poisson method). The corrected PDT is for both methods, the PDT minus the average of all background data (<PDT).

The calculation of the limit of detection (LOD) for particle masses is presented in eqn (8) (adapted from ref. 32),

and the calculation of the limit of detection (LOD) for the spherical particle sizes as expected size equivalents in eqn (9) (adapted from ref. 32),

Gaussian peak-fitting. A common issue of spICP-MS is the accurate positioning of the PDT. Even if the PDT is set well, some background events may remain and would be combined with the particle events, leading to an underestimation of the particle sizes when an averaging of the overall dataset is used. Thus, Sparta uses a flexible multi-modal Gaussian peak-fitting algorithm (based on ref. 45) to simplify the identification of detected significant masses or sizes, which can account for one to four Gaussian peaks which are either particle distributions or the background itself (e.g., Fig. 3a). In addition, the alpha (α) and beta (β) errors of the identified peaks are included in the code. The α error provides the confidence interval of the calculated peak values (height and deviation). The β error provides the probability of a Type II error (β) assuming the true amplitude is the estimated one. In our case, we obtain an indication of the peak detectability based on the significance of the amplitude (i.e., using the fitted amplitude as a proxy for the true effect size of the β values), providing in return the power of the test (1−β). In this way, both parameters (α and the power) allow proper identification of Gaussian peaks with a default value of 95% (k = 2) confidence (i.e., modifiable to higher reliabilities with the ‘alpha_error’ parameter in the code). All masses and sizes were plotted as histograms, and after peak-fitting it was statistically decided (via the power) which peaks are real and suitable (e.g., case of peak-fitting for SiO₂, 300 nm, Agilent, 0.1 ms dwell time in Fig. 3a and b). In several cases, peak-fitting improved the accuracy of the mass and size calculations (e.g., Fig. 3b, orange vs. blue boxplot), otherwise underestimated due to spurious NPs above PDT emerging from the background. The Gaussian peak-fitting algorithm extracts the original data (blue bars) under each detected Gaussian peak (orange), even if the bar frequency (no. of events) exceeds the orange fit.

Fig. 3 Gaussian peak-fitting algorithm for silica (SiO₂), 300 nm, 0.1 ms dwell time (Agilent). (a) Histogram of all particle events (blue bars) above particle detection threshold (PDT) and Gaussian peak-fitting (orange, n = 3) for the Poisson method. Peaks marked in ‘bold’ (i.e., peak 1) represent ‘real’ peaks, which are identified via the power test of the β error, meaning they are within the confidence interval of 95% (k = 2). Other peaks (i.e., peak 0 and 2) are below the confidence interval of 95% and thus marked in brackets and not bold. The red dashed line represents the LOD_size. The standard deviation (SD) emerges from the Gaussian function (k = 1). (b) The blue boxplot shows all particle events above PDT and the orange boxplot shows only the extracted particle events under the Gaussian peak-fitting. The certified diameter (average, green line ± uncertainty) was analysed using transmission electron microscopy (TEM) with a certified uncertainty (k = 1).

Results and discussion

Size method comparison and validation

In this section, we validate Sparta in a benchmarking round using NP sizes of three different elements/species (Au, TiO₂ and SiO₂) measured with two instruments (Agilent and Thermo). As the NP sizes emerge from the NP masses, we only present the obtained sizes. Examples of size histograms from Sparta can be found in Fig. S2–S6. In the following sections, the terms ‘Agilent’ and ‘Thermo’ are used for the instruments 8900 ICP-MS/MS and Thermo iCAP-TQ, as specified in section ‘spICP-MS instrumentation’.

Gold (Au) – ∼30 nm. One of the most studied NPs using spICP-MS are Au-containing NPs, as Au has a low (ionic) background and little interferences. Au is also often used to determine the transport efficiency in spICP-MS. The manufacturer software (MassHunter from Agilent and Qtegra from Thermo) usually shows excellent results for Au, this is why we used Au as a validation tool for Sparta. As hypothesised, all methods, dwell times and instruments showed very similar results close to the certified nominal diameter of 30 nm (black dashed horizontal line) (Fig. 4). Only the 0.1 ms results fit better to the certified mode of 32.7 nm (green horizontal line, Fig. 4), corresponding to the hydrodynamic diameter certified obtained via nanoparticle tracking analysis (NTA) from the certificate.

Fig. 4 Size comparison for the ∼30 nm Au NPs via spICP-MS between Agilent (0.1 vs. 5 ms dwell time) and Thermo (5 ms dwell time) using different data processing methods. The indicative hydrodynamic diameter (mode, green line) ± uncertainty (green area) is taken from the providers' certificate, originally measured via nanoparticle tracking analysis (NTA) with an expanded uncertainty (k = 2). The black dashed line represents the nominal particle diameter of 30 nm from the provider. NP-events were identified above the particle detection threshold (PDT), defined as follows: for the MassHunter software (Agilent), options were set on ‘auto” (algorithm not known), for SPCal (set on ‘auto’, based on the background) the chosen method was Poisson (both instruments and dwell times), and for the Excel (both instruments), the Gaussian method (average (µ) + 6 standard deviations (SDs)) was used. In our case, Sparta used the Poisson method (all data above PDT) for both instruments and dwell times. The Qtegra (Thermo) software was not applied in this round of measurements. ‘IQR’ = interquartile range. ‘LOD’ = size limit of detection. ‘BED’ = background equivalent diameter.

For Sparta, we used the Poisson method to set the particle detection threshold (PDT) and identify NP events. Anyhow, there were only minor differences between both methods, as an example for the case of 0.1 ms dwell time via Agilent, the median for the Gaussian method (µ + 7SD) was 30.9 nm vs. 31.2 nm for the Poisson method. The Gaussian peak-fitting is not shown here, instead the NP events represent all valid events (no agglomerates) above PDT. The size limit of detection (LOD_size) and the background equivalent diameter (BED) are lower for a short dwell time of 0.1 ms compared to 5 ms, suggesting improved particle discrimination from the background for shorter dwell times. Sparta shows a LOD_size of ∼6.6 nm for Agilent at 0.1 ms dwell time, whereas at 5 ms it increased to almost double (∼12.9 nm) for Agilent and nearly triple (∼16.6 nm) for Thermo. This is explained by the fact that with increasing dwell time, a higher background signal is integrated within a specific dwell time, resulting in a lower signal-to-noise ratio.²⁴ Lower LOD_size with decreasing dwell time is a common feature reported in the literature, e.g., LOD_size = 9.0 nm (0.1 ms) vs. 12.7 nm (5 ms) for Au NPs.^2,3 For 0.1 ms dwell time, our LOD_size (via Sparta) is very similar to the one presented in the Agilent application note (∼6.5 nm),²⁹ and the LOD_size for 5 ms dwell time is for Agilent very close and for Thermo slightly larger than reported in literature.^2,3 In a nutshell, for 30 nm Au NPs, the applied instruments, data processing methods and dwell times only show minor differences in terms of their calculated particle size and work well.

Titanium dioxide (TiO₂) – ∼74 nm. TiO₂ NPs are more challenging than Au for spICP-MS. Although reported LOD_size (0.1 ms dwell time) is ∼30 nm, i.e., relatively higher than Au NPs,^2,30 our LOD_size for 0.1 ms dwell time was ∼43 nm with Sparta using the Poisson method (Fig. 5). The resulting median size obtained via Sparta (∼71 nm, 0.1 ms dwell time) was close to the certified size (74 ± 10 nm). For the same dwell time, both MassHunter (median ∼96 nm) and SPCal (median ∼209 nm) overestimated the certified size. For SPCal, it seems that the threshold was set too high so that the small events were lost within the background. This was supported by the low number of events detected (n = 29) and the comparably high LOD_size (∼63 nm).

Fig. 5 Size comparison for the ∼74 nm TiO₂ NPs via spICP-MS between Agilent (0.1 vs. 5 ms dwell time) and Thermo (5 ms dwell time) using different data processing methods. Two diameters (average, green line) ± uncertainty (green area) were used from the provider's certificate obtained using laser diffraction spectrometry (71 ± 4 nm) and X-ray disc centrifugation (77 ± 7 nm) with a certified expanded uncertainty (k = 2). NP-events were identified above the particle detection threshold (PDT), defined as follows: for the MassHunter software (Agilent), options were set on ‘auto’ (algorithm not known), for SPCal (set on ‘auto’, based on the background) the chosen method was Poisson (both instruments and dwell times), and for the Excel (both instruments), the Gaussian method (average (µ) + 7 standard deviations (SDs)) was used. In our case, Sparta used the Poisson method (all data above PDT) for both instruments and dwell times. The Qtegra (Thermo) software was not applied in this round of measurements. ‘IQR’ = interquartile range. ‘LOD’ = size limit of detection. ‘BED’ = background equivalent diameter.

For 5 ms dwell time, previous literature reported a LOD_size of ca. 42 nm (measured as ⁴⁷Ti).² However, we used the isotope ⁴⁸Ti as we do not expect to have the ⁴⁸Ca interference and ⁴⁸Ti has a higher abundance, leading to a higher sensitivity. All methods with the 5 ms dwell time show very similar results, systematically overestimating (median: ∼90–110 nm) the certified particle size. This suggests that smaller NPs were lost within the background for 5 ms dwell time, which is an instrument, method-based issue, not a data processing effect. This is further confirmed by the trend of a higher LOD_size with increasing dwell time, already seen for Au NPs (Fig. 4). This is clear from Sparta, where the LOD_size is higher than the certified NP size (Sparta, 5 ms dwell time: LOD_size ∼82 nm). The lower LOD_size for SPCal and lower BED for MassHunter is a result of different calculations. SPCal uses their PDT to obtain their LOD_size,³² instead of our corrected PDT (both are different), as explained in Subsection ‘Particle detection threshold (PDT)’ in the Materials and methods. In contrast, MassHunter has a different approach to estimate their BED, based on an average of their particle signals and the intensity of noise signals, which is defined in detail in ref. 46. Moreover, minor size differences were seen between the Thermo and Agilent instruments for 5 ms dwell time. In summary, when measuring TiO₂ NPs ≤ 82 nm, a shorter dwell time, such as 0.1 ms, should be preferred to accurately determine the NP sizes. Under these conditions, Sparta performs well, i.e., was successfully validated using TiO₂ NPs for 0.1 ms dwell time.

Silica (SiO₂) – ∼50 nm. Even more challenging than TiO₂ are SiO₂ NPs for spICP-MS, especially small sizes. As mentioned in the introduction, SiO₂ suffers from polyatomic interferences and Si-based NPs from glass materials equipped in the instrument, resulting in comparable poor signal-to-noise ratios.^13,17 In this work, the smallest particle size studied was ca. 50 nm SiO₂, certified by the provider via TEM (Fig. 6). Results show that for the short dwell time (0.1 ms), all methods applied (MassHunter, SPCal and Sparta) systematically overestimate (<100 nm) the certified diameter. Nevertheless, out of all methods, Sparta shows the closest size compared to the certified diameter (median ∼65.3 nm). The explanation for the overall overestimation is that small NPs cannot be differentiated from the background and are thus not identified as such. An earlier publication also found that 80 nm SiO₂ NPs were partly lost within the background and the size was therefore overestimated.¹⁷ To compensate for the overestimation, the authors used a Gaussian function, overlapping with the background, to accurately determine the particle distribution. Therefore, small SiO₂ NPs seem to be still a limitation for the spICP-MS technique, even if choosing small dwell times (i.e., 0.1 ms).

Fig. 6 Size comparison for the ∼50 nm SiO₂ NPs via spICP-MS between Agilent (0.1 vs. 5 ms dwell time) and Thermo (5 ms dwell time) using different data processing methods. The certified diameter (average, green line ± uncertainty) is given by the provider via transmission electron microscopy (TEM) with a certified uncertainty (k = 1). NP-events were identified above the particle detection threshold (PDT), defined as follows: for the MassHunter software (Agilent), options were set on ‘auto’ (algorithm not known), for SPCal (set on ‘auto’, based on the background) the chosen method was for 0.1 ms Gaussian (average (µ) + 5 standard deviations (SDs)) but changed to Poisson as this was the more suitable approach here and for 5 ms dwell time (both instruments) the Gaussian (µ + 5SD) was chosen. For the Excel (both instruments), 5 ms dwell time, the Gaussian method (µ + 4SD) was used. In our case, Sparta used the Poisson method (all data above PDT) for 0.1 ms and for 5 ms dwell time, the Gaussian method (µ + 5SD, all data above PDT) was used for both instruments. The Qtegra (Thermo) software was not applied in this round of measurements. ‘IQR’ = interquartile range. ‘LOD’ = size limit of detection. ‘BED’ = background equivalent diameter.

In contrast, the results achieved using a higher dwell time of 5 ms greatly overestimate the particle sizes with all data processing methods and both instruments. This becomes even clearer when emphasising the lower total number of events detected at 5 ms dwell time compared to 0.1 ms dwell time, although the total acquisition time at 5 ms was three times larger than at 0.1 ms dwell time (120 s vs. 40 s). In any case, a high number of events is sought to achieve sufficient statistics. A minimum number of approximately 100 events per minute acquisition time is suggested but may be set lower for homogeneous materials.⁴⁷ In the case of 0.1 ms dwell time (40 s acquisition time), ∼66 NPs are suggested, which is reached by Sparta (282 events). Under these conditions, the MassHunter software (Agilent) shows the closest results to the certified size for a 5 ms dwell time. However, it is highly unlikely to be true, as the BED was higher than the greatest event, i.e., the largest NP, detected. Given the fact that the exact MassHunter software algorithm to determine the PDT is not disclosed, we cannot verify this hypothesis, nor explain the reason for the good output for a 5 ms dwell time. In any case, even if the diameter for 50 nm SiO₂ NPs is overestimated for both dwell times, Sparta shows comparable, or even better, results to other data processing methods, representing the best-fitting diameter compared to the certified value for 0.1 ms dwell time.

Silica (SiO₂) – ∼100 nm. By increasing the NP diameter by a factor of two, for 0.1 ms dwell time, all three data processing methods show very well-fitting results compared to the certified diameter (Fig. 7). This is a big step forward compared to previous literature, showing that only 1–3% of 100 nm SiO₂ NPs could be detected, although this was using dwell times of 25–100 µs, and explained that only the largest particles were differentiated from the background.⁸ However, as with the 50 nm SiO₂ NPs, the results for 5 ms dwell time still overestimate the certified diameter. MassHunter (Agilent) gives once more a well-suited diameter, even though the BED is higher than the greatest event detected. However, for the first time, we detect differences between the two instruments (Agilent vs. Thermo). The data processing of the output from the Agilent instrument results in many NPs events with the SPCal, Excel and Sparta methods, despite the difference between 0.1 ms and 5 ms analyses (i.e., fewer events in 5 ms), indicating that only larger particles were detected using 5 ms dwell times and small NPs were lost in the background. In contrast, the number of detected events obtained from the Thermo instrument was much lower, suggesting that even more (small) NP events were lost in the background. Nevertheless, despite these NP losses, the overall NP sizes calculated from 5 ms dwell times are systematically overestimated with both instruments and generally with all data processing methods. Moreover, in this case, the Thermo software (Qtegra) displayed an average ± SD (219 ± 14 nm), which also overestimates the size, indicating that the calculations of the commercial software are similar to SPCal, Excel and Sparta. Thus, as expected, the dwell time of 0.1 ms gave the most accurate NP size calculations for SiO₂ ∼100 nm, almost independent of the data processing method.

Fig. 7 Size comparison for the ∼100 nm SiO₂ NPs via spICP-MS between Agilent (0.1 vs. 5 ms dwell time) and Thermo (5 ms dwell time) using different data processing methods. The certified diameter (average, green line ± uncertainty) was analysed using transmission electron microscopy (TEM) with a certified uncertainty (k = 1). NP-events were identified above the particle detection threshold (PDT), defined as follows: for the MassHunter software (Agilent) and Qtegra software (Thermo), options were set on ‘auto’ (algorithm not known), for SPCal (set on ‘auto’, based on the background) the chosen method was Gaussian (average (µ) + 5 standard deviations (SDs)) for both dwell times and instruments, and for the Excel (both instruments), the Gaussian method (µ + 4SD) was used for 5 ms dwell time. In our case, Sparta, the Gaussian method (µ + 7SD, all data above PDT) was used for 0.1 ms and the Gaussian method (µ + 5SD, all data above PDT) for 5 ms dwell time (both instruments). ‘IQR’ = interquartile range. ‘LOD’ = size limit of detection. ‘BED’ = background equivalent diameter.

Silica (SiO₂) – ∼300 nm. Out of all tested particle types and sizes, the outcomes from the measurements of SiO₂ NPs of ∼300 nm showed best the advantages of the newly developed Gaussian peak-fitting/extraction algorithm integrated into Sparta, applied here for both dwell times and instruments. For the short dwell time (0.1 ms), we obtained a median ± SD of 279.1 ± 17.1 nm (Fig. 8). Both, MassHunter and SPCal calculations underestimated the certified particle diameter by a factor of 2–3, the same factor as Sparta if the peak-fitting/extraction algorithm, as presented in Fig. 3, would not have been applied. As discussed, for the ∼50 nm SiO₂ NPs, a higher event number (n) usually provides better statistics. However, this does not necessarily mean that a higher event number also provides a better data quality. As demonstrated in Fig. 3, all events above PDT would result in many false-positive background signals (n = 2091) detected as spurious events, vs. the Gaussian peak-fitting algorithm provided fewer events (n = 255) but a more accurate size. The higher background (originally in ethanol matrix) in the case of 300 nm SiO₂ might be related to a lower dilution factor (DF: ∼6 × 10⁶) as the particle concentration was lower compared to 50 nm (DF ∼1 × 10⁹) and 100 nm SiO₂ (DF ∼2 × 10⁸). As the original silica standards were provided from the manufacturer in glass vessels, which were bath sonicated to deagglomerate NPs before dilution, there are most likely more silica NPs emerging from the vessel itself. A lower dilution factor directly means more silica NPs in the final sample.

Fig. 8 Size comparison for the ∼300 nm SiO₂ NPs via spICP-MS between Agilent (0.1 vs. 5 ms dwell time) and Thermo (5 ms dwell time) using different data processing methods. The certified diameter (average, green line ± uncertainty) was analysed using transmission electron microscopy (TEM) with a certified uncertainty (k = 1). NP-events were identified above the particle detection threshold (PDT), defined as follows: for the MassHunter software (Agilent) and Qtegra software (Thermo), options were set on ‘auto’ (algorithm not known), for SPCal (set on ‘auto’, based on the background) the chosen method was Gaussian (average (µ) + 5 standard deviations (SDs)) for both dwell times and instruments. For Excel, the Gaussian method (µ + 6SD for Agilent and µ + 4SD for Thermo) was used for a 5 ms dwell time. In our case, Sparta, the Poisson method (peak-fitting/extraction algorithm) was used for 0.1 ms and the Gaussian method (µ + 6SD for Agilent vs. µ + 5SD for Thermo) was used for 5 ms dwell time, both using the peak-fitting/extraction algorithm. ‘IQR’ = interquartile range. ‘LOD’ = size limit of detection. ‘BED’ = background equivalent diameter.

It is worth noting that 5 ms dwell time measurements only show acceptable sizes for larger SiO₂ NPs, well-fitting for both instruments when using SPCal, Excel and Sparta as data processing methods. In contrast to Sparta, in combination with the peak-fitting/extraction algorithm, Excel takes all data above PDT. For a 5 ms dwell time, the LOD_size was much higher (LOD_{size, Sparta}: 203 nm) than for 0.1 ms dwell time (LOD_{size, Sparta}: 60 nm). Thus, the small spurious background particles, which were only detected at 0.1 ms dwell time, played only a minor role at 5 ms dwell time. Nevertheless, even at 5 ms dwell time, Sparta using the peak-fitting/extraction algorithm, was able to remove some leftover background, which resulted in better fitting sizes (median ± SD: 279.4 ± 20.1 nm). Both manufacturer software (MassHunter and Qtegra) slightly underestimated the particle diameter. Overall, we show that Sparta also works well for larger SiO₂ NPs (∼300 nm) with both dwell times and instruments and highlight the importance of the Gaussian peak-fitting/extraction algorithm.

Particle number concentration method comparison and validation

A second important parameter, which was used alongside the sizes for the benchmarking and validation of Sparta, is the element-specific particle number concentration (PNC), which shows the uniqueness of spICP-MS compared to other methods such as Nanoparticle Tracking Analysis or Laser Induced Breakdown Detection. In this study, we aimed to target an accuracy of one order of magnitude in PNC. For Au, the outcome is similar as for the sizes, the PNC can be accurately estimated with any instrument, dwell time and data processing method, including Sparta, which suits well to the other methods and was thus successfully validated (Fig. 9a). The error is greater for TiO₂, even when using the more sensitive 0.1 ms dwell time, however, all methods, instruments and dwell times were within one order of magnitude of the calculated PNC assuming 74 nm as an average diameter, as the PNC was not certified here (eqn (S3) and (S4)). Thus, Sparta was also successfully validated for TiO₂.

Fig. 9 Particle number concentration comparison (semi-logarithmic) via spICP-MS between Agilent (0.1 vs. 5 ms dwell time) and Thermo (5 ms dwell time) using different data processing methods. (a) Au, ∼30 nm and TiO₂ ∼74 nm (differences in the calculated PNCs between Agilent and Thermo are a result of two different PNCs of the original suspensions). (b) SiO₂, ∼50, ∼100 and ∼300 nm. For spICP-MS, data represent the average ± 2SD of three replicates. For ‘SiO₂, 50 and 100 nm via Thermo, 5 ms (all methods)’ and ‘SiO₂, 50 nm via Agilent (MassHunter), 5 ms’, only two replicates were available, thus the error represents half the range. The certified PNC (average, green line ± uncertainty, k = 2) was certified using nanoparticle tracking analysis (NTA). NP-events were identified above the particle detection threshold (PDT), which were further specified in the descriptions of Fig. 4–8.

As expected, the accurate determination of the PNC of SiO₂ NPs (Fig. 9b) was more complex and challenging than for Au and TiO₂. When comparing the results from the short dwell time (0.1 ms), the certified value of the PNC for the 50 nm NPs is slightly underestimated, while for 100 and 300 nm it is rather overestimated. This may be because, for 50 nm, very small particles may be lost within the background, while for 100 and 300 nm, the detected background (i.e., polyatomic interferences, glass containing NPs) increases the actual number of detected events, thus increasing the total PNC. This issue was resolved for 300 nm SiO₂ and 0.1 ms, because the peak-fitting/extraction algorithm was used to remove spurious NP events, leading to a well-fitting PNC. Nevertheless, all three data processing methods are within the accuracy of one order of magnitude in PNC.

In contrast, the longer dwell time (5 ms) shows an underestimation of PNC for the 50 nm NPs with all methods for Agilent and even more for Thermo. This can be explained, on the one hand, by a lower LOD_size achieved by Agilent (more NPs detected) compared to the Thermo instrument. On the other hand, for the Thermo instrument, a cyclonic spray chamber with a lower η was used (∼5%), which has almost no effect on the NP mass or size calculation for homogeneous samples but the uncertainties in η lead to a larger relative error in PNC, especially for smaller η as is here the case. The underestimation in PNC achieved by Thermo is also visible for the 100 nm NPs, in contrast to the well-fitting results achieved by the Agilent instrument via SPCal, Excel and Sparta, despite the worse fitting for its own software (MassHunter). We believe that for Thermo, the LOD_size was too high to accurately estimate the PNC here, receiving underestimated results with all methods. For the 300 nm SiO₂ NPs, Sparta suits best to the certified PNC, together with the Thermo software (Qtegra). The difference between the instruments does not play a big role here as the LOD_size is, in this case, lower than the actual NP sizes, even though Excel still significantly underestimates the PNC for the Thermo instrument.

In summary, we successfully validated Sparta also for SiO₂ in PNC and could accurately estimate the PNC within one order of magnitude for all measurements with 0.1 ms dwell time and for 5 ms dwell time for Agilent ≥ 100 nm and Thermo ≥ 300 nm.

Working features, new tool improvements and remaining limitations of Sparta

Our proposed Python data processing algorithm (Sparta) was applied to model NPs above and compared and validated using two different instruments and four (non-)commercial data processing software. This section contains a general summary and discussion of the applied working features, the new improved tools and the remaining limitations which should be tackled in future studies. An overview is presented in Table 2:

Table 2 Overview of the new features/tools of Sparta, their improvement potential and original sources

Feature/tool	New	Improvement	Original source(s)
Baseline drift correction	Calculation (mode instead of average; eqn (1))	More robustness using mode instead of averages to directly access background and not particles. No need for dynamic PDTs if noise is consistent (baseline drift but no noise width fluctuations), only one PDT for the whole dataset	34
Particle detection threshold (Gaussian method; µ + 3–7 SD)	(1) ‘clear_cut-function’. (2) Buffer of 40 iterations	(1) Optional feature: no background contribution for monodispersed well-defined NPs (e.g. Au). (2) Ensured stable PDT	6 and 37
Particle detection threshold & critical value (Poisson method)	Adding the average from all raw data (eqn (2) and (3))	Allowing to use the Poisson method also for high background systems (e.g. silica)	32, 38, 42 and 43
Peak summation	Peak summation for datasets if choosing dwell times <2 ms	Necessary, as each NP event may have several data points	32
Outlier removal	Outlier removal from valid NP events (>µ + i SD)	Optional feature: Elimination of false-positive artefacts or agglomerates for well-defined, monodisperse CRMs	—
Transport efficiency (particle number/frequency method)	Calculation slightly adapted; (eqn (4))	—	6, 32 and 41
Particle number concentration	Possibility to calculate PNC from an extracted peak using the Gaussian peak-fitting/extraction algorithm (eqn (5))	PNC calculation of artefact-free and/or polydisperse distributions	32
Particle masses/sizes	Calculations slightly adapted (eqn (6) and (7))	(1) Mass fraction moved to the particle size calculation to report particle masses for unknown materials without any assumptions. (2) consideration of the intercept of the ionic, linear calibration	32
LOD	Corrected PDT (PDT minus the average of background data); (eqn (8) and (9))	Consideration of the intercept from the ionic, linear calibration	32
Gaussian peak-fitting/extraction	Multi-modal peak finding, fitting and extraction algorithm for 1–4 peaks within the NP histogram for NP masses and sizes. Including alpha (α)/beta (β) error calculations for statistical significance of the identified peaks	Removal of remaining, false-positive background events; investigation of chosen peaks from polydisperse distributions. Peak descriptors are better constrained with the α error calculation. Additionally, the power test of the β error calculation improves the confidence of peak identifications	—

---

Baseline drift correction. The potential baseline drift was corrected as a very first step in Sparta before further data processing. As presented in Fig. 1, our approach allows the usage of only one PDT for the whole dataset. Previous literature³² described local data thresholding to correct for baseline drift scenarios, delivering several dynamic PDTs, which worked well but could lead to confusion due to several PDTs during data interpretation. If not correcting baseline drift and also not considering dynamic PDTs, the particle mass/size distributions but also PNCs are directly affected. In the case of a negative baseline drift (e.g. Fig. 1a, ∼15–30 seconds), particle events would be lost within the background and detected events would be underestimated due to a lower signal intensity. In the opposite case, having a positive baseline drift (e.g. Fig. 1a, ∼5–15 seconds), false-positive background events would be detected as particles and particle mass/size distributions would be overestimated due to a higher signal intensity. Compared to the original idea introduced by Cornelis & Rauch 2016, we used moving modes instead of moving averages, which are more robust against extreme values (i.e., particle events) to directly access the background signals and their drift.³⁴ Our baseline drift correction was developed and validated for rather small baseline fluctuations such as pumping inconsistencies within one sample, but not for changes in the noise width (e.g., ref. 34). Overall, we did not discover large changes in the noise width in all our samples and the baseline drift correction homogenises the time-dependent variations of the spectrogram, keeping the proportionality between the ionic signal and the detectable peak events.

Particle detection threshold. For the Gaussian method, we first developed an optional feature applicable for well-defined, monodisperse CRMs such as Au NPs, which detects a so-called ‘clear_cut’ between the background and the NP distribution within the histogram. It is an automatic feature without the user-dependent ‘error contribution’. If enabled, it would replace the statistical Gaussian threshold, which suffers in some cases from remaining, false-positive background events. Secondly, we implemented a buffer of 40 iterations, as more iterations do not change the final PDT and do not take up much time, although in most cases, fewer than ten iterations were sufficient to define the Gaussian PDT. For the Poisson method, we adjusted the original Poisson formulas (critical value and PDT) as described in Currie 1968, by adding the average of the raw data.³⁸ The original Poisson approach is specifically applied and validated for low-intensity systems (e.g. low background such as for Au).³³ This approach increases the PDT for low and medium background systems such as Au and TiO₂, respectively, consequently also the LOD increases. In any case, the effect is minor, as the LOD only increases for Au from 6.4 to 6.6 nm and for TiO₂ from 42.0 to 43.4 nm (for Agilent, 0.1 ms dwell time). However, our new Poisson approach now allows data processing also for high background systems such as SiO₂ (as presented in Fig. 6 for 0.1 ms dwell time, Agilent). In contrast, the conventional Poisson threshold (Fig. S7a) would have detected several false positive background signals and was not suitable, while the adapted approach from Sparta (Fig. S7b) shows well-fitting PDT_Poisson. In summary, our empirically tested adjustment accounts for a more realistic detection scenario, where the threshold must exceed both background average and the statistical variation (noise). The workability was presented and validated in the above sections, comparing our results with other the comparison methods.

Peak summation. The peak summation (only activated if the dwell time is <2 ms) works that way that adjacent data point (here referred as peak) above the background threshold (i.e., the PDT value) are added continuously until the upcoming peak trend (i.e., the following data point) falls below the PDT (Gaussian method) or below the critical value (Poisson method). Like the area under a curve, approximated as the sum of the individual bar charts. Each individual peak height is registered without the contribution of the ionic background, as the ionic average is subtracted from the individual peak heights before adding them up. Finally, the summed peak (in ‘counts per seconds’) represents one detected particle event. This operation is crucial, as otherwise the particle number concentration would be overestimated, and their mass and size would be underestimated, as ‘combined particle events’ would have been split into individual events. In contrast, when the dwell time is >2 ms, adjacent peaks (above the PDT) are considered as independent particle events (i.e., 5 ms is in most of the cases larger than the event duration), based on good sample dilution and thus, the peaks are accounted for as individual peaks for the total peak budget. Moreover, in Fig. 3–8, the variable ‘n’ represents the number of events detected by each data processing method. As each method works slightly differently (e.g. in determining the PDT and executing the peak summation), there are differences for ‘n’ visible even within the same dataset and dwell time. Comparing both dwell times, lower dwell time means we are able to detect smaller NPs and therefore resulting in more events if the LOD of 5 ms dwell time is close or higher than the certified size (e.g. 50 and 100 nm SiO₂, Fig. 6 and 7). However, the differences in total acquisition time must be considered if directly comparing event numbers (for 5 ms was 120 s and for 0.1 ms 40 s). It should be further noted that the peak summation was only applied for the short dwell time, necessary if the dwell time is shorter than the NP duration. Finally, event number differences between the instruments can occur due to daily varying instrument sensitivity.

Outlier removal. This tool was used for our well-defined, monodisperse model NPs to remove outliers (artefacts/agglomerates) from valid particle events above PDT, which would potentially overestimate the particle masses and sizes. For the PNC calculation, all events (also outliers) are considered as they represent still particle events, although if agglomerated, the PNC would be slightly underestimated. If particles are normally distributed, no outliers would be detected and none removed. In any case, detected outliers were very few, not reproducible and in most cases <5 events (e.g. Fig. S1). Nevertheless, we recommend disabling this tool (Sparta version 1.1, disabled by default) for natural, unknown and/or heterogeneous materials, as the tool would then remove actual large particle events.

Transport efficiency. In this study, we used and evaluated the particle number/frequency method, as the PNC was certified for the 30 nm Au NPs, but the size was determined using NTA, representing the hydrodynamic diameter. Therefore, it was more accurate to use the proposed method. Generally, also the particle size, waste collection or dynamic mass flow methods can be used to determine η.^23,41,48 Briefly, the particle number method counts the events in a given time and compares them with the total number of events (known/certified) in the suspension by using the volume flow for the calculation.⁴⁸ In contrast, the particle size method uses the slope of the ionic calibration, the volume flow and dwell time to convert the mean particle intensity to an elemental mass and compares it to a certified elemental particle mass.⁴⁸ Both methods require CRMs (e.g. Au NPs).⁴⁸ The waste collection method determines gravimetrically the sample inlet and waste outlet flow to estimate the ratio of suspension arriving in the plasma, but could not be applied here as we used a total consumption spray chamber.⁴⁸ The waste collection method has further the disadvantage that adsorbed droplets in the tubes or on the spray chamber, as well as evaporation, are not considered, which leads to a large uncertainty that is ∼20% relative standard deviation (RSD).⁴⁸ Lastly, the dynamic mass flow method also works gravimetrically and compares the slope of the direct continuous sample mass flow measurements over time to the slope of the uptake by ICP-MS.^23,48,49 The latter two methods have the advantage that they do not require CRMs. In any case, most authors recommend either the particle size or particle number method, depending on each special case.^23,48 The suitability of the particle number/frequency method, compared to the particle size and waste collection method, was reported in a previous publication.¹³

Particle number concentration. The PNC was discussed in detail in the above section. In Sparta, we added the possibility to calculate the PNC within chosen detected peaks to accurately assess PNCs in polydisperse distributions. This was shown and validated for the case of 300 nm SiO₂ NPs (0.1 ms), where the PNC was determined via the Gaussian peak-fitting algorithm using the extracted events under the detected peak, without the remaining background events (Fig. 9b).

Particle masses, sizes and their LOD. The particle sizes emerge from the particle masses and were discussed for each material accordingly. The only major adaptation made in (eqn (6)–(9)) is the consideration of the intercept of the linear, ionic calibration, as we find this is part of the equation and should be considered. Furthermore, we moved the mass fraction to the particle size calculation, as for e.g. natural/unknown materials, no mass fraction is usually known and the particle masses can then still be used without any assumptions. For the LOD_mass, we used the corrected PDT, which means the original PDT minus the average of background data (all data below PDT), as the same operation was also done for the raw particle intensities. The LOD_size emerges from the LOD_mass.

Gaussian peak-fitting/extraction. This tool was applied for the case of 300 nm SiO₂ NPs (0.1 ms) and discussed in detail in the above section. This is especially valuable to remove false-positive background contributions or to determine specific peaks within polydisperse distributions regarding their mass, size and/or PNCs. This approach is reinforced by the statistical calculation of the α and β errors of the peak identifications, providing a higher confidence (here 95%, k = 2) on the NP identification.

Conclusions

In this study, we present a novel and transparent data processing Python algorithm (Sparta) capable of processing quadrupole-based single particle Inductively Coupled Plasma-Mass Spectrometry (spICP-MS) raw data from different manufacturers (here Agilent Technologies and Thermo Fisher Scientific). The elaborated benchmarking of two instruments, micro- and millisecond dwell times as well as the use of five data processing methods not only validates Sparta, but also shows the improvement in accuracy, especially reached for challenging silica (SiO₂) nanoparticles (NPs) in size and particle number concentration (PNC).

We present that all data processing methods with both instruments provide well-fitting sizes and PNCs for 30 nm Au NPs, further validating Sparta. Sparta is a step-forward in data processing as it combines previous well-working procedures with new ideas and adaptations, especially well presented for microsecond dwell times, showing accurate sizes and PNCs up to an order of magnitude even for matrix-challenging materials, i.e., TiO₂ and SiO₂ NPs. In general, shorter dwell times (0.1 ms) significantly improve the limit of detection for sizes (LOD_size), allowing the accurate detection of 74 nm TiO₂ as well as 100 and 300 nm SiO₂ NPs, although 50 nm SiO₂ NPs were too close to the LOD_size and therefore overestimated with all data processing methods and instruments. One highlight is the baseline-correction, which allows the usage of only one PDT for the whole dataset and thus solves errors when subtracting the average background signal from each single particle event. The microsecond dwell time approach is combined with a peak summation to correctly identify separate particle events, with the advantage of decreasing the signal-to-noise ratio and the size limit of detection of NPs compared to millisecond dwell times. Furthermore, a widely discussed issue is the threshold method and criterion used for the correct positioning of the particle detection threshold (PDT). Usually, setting the PDT is a compromise between the reduction of false-positive NP events originating from the background (higher PDT) vs. the prevention of lost events in the background (lower PDT). Sparta provides the flexibility of using either the Poisson method (usually for low background elements such as Au) or choosing the conventional Gaussian method (usually for high background elements such as Si) to determine the PDT. Events above the PDT can be further divided into peaks using the Gaussian peak-fitting/extraction algorithm to eliminate spurious ‘leftover’ background signals, as presented in this study for the 300 nm SiO₂ NPs. Consequently, a too-low PDT is no longer problematic for well-known materials when applying the Gaussian peak-fitting/extraction algorithm, as false-positive NP events are filtered out after setting the PDT. Moreover, the Gaussian peak-fitting allows the extraction of chosen size distributions for polydisperse samples, which was previously not possible as all events above PDT were extracted and merged to the result. For unknown heterogeneous samples, the correct setting of the PDT remains important, particularly when analysing small NPs whose distribution might overlap with the background.

Beyond that, Sparta can not only be used for technical, engineered NPs, but can be readily applied to all other NPs of technical and natural origin. For instance, the quantification of the element association of metals into/onto single biological cells (single cell ICP-MS) or even the elemental adsorption onto nano or microplastics can be studied with Sparta (already tested, not shown here). Combined with the Gaussian peak-fitting/extraction algorithm, microplastics can be even separately detected and differentiated from small organic colloids when tracking carbon. Both applications are becoming increasingly important to better understand uptake and/or transport mechanisms of pollutants, e.g., in aquatic environments, thus validated analyses and data processing methods, as the one shown in this work, are of utmost relevance to facilitate future studies and their comparability in this field.

Author contributions

Conceptualisation, T. S., M. M.-B.; methodology, S. H.; validation, S. H., T. G. D., M. C.-R.; investigation, S. H., M. C.-R.; coding, S. H. and T. G. D.; resources, T. S., M. M.-B.; writing – original draft preparation, S. H.; writing – review and editing, S. H., T. G. D., D. M., M. M.-B. and T. S.; visualisation, S. H.; supervision, D. M., T. G. D, M. M.-B. and T. S.; all authors have read and agreed to the published version of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Data availability

Code availability: The code is available via GitHub for reference materials (e.g. Au) to calculate the transport efficiency as well as a code for all other elements as an example for SiO₂ under the following link: https://github.com/Steffen9204/Sparta-sp_scICP-MS. The version of the code employed for this study is Version 1.1.

The datasets used can be found in the supplementary information (SI) or are available from the corresponding author upon reasonable request. Supplementary information is available. See DOI: https://doi.org/10.1039/d5ja00285k.

Acknowledgements

S. H. would like to acknowledge the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) (Project SCHA 1854/7-1, DFG SPP2005, Opus Fluidum Futurum and INST 275/441-1 FUGG, ICP-MS/MS) for funding this research. S. H. would further like to thank David Clases for taking the time to give insights into spICP-MS data processing during the ‘ICP-Anwender*innen Treffen’ conference (2022) in Leoben (Austria), and for providing and further developing the freely available software SPCal, which helped us to advance Sparta. Finally, S. H. appreciates the contribution of David Ho, who helped with the initial stages of the Python-algorithm creation.

References

1. C. Degueldre and P. Y. Favarger, Colloids Surf., A, 2003, 217, 137–142.
2. F. Laborda, A. C. Gimenez-Ingalaturre, E. Bolea and J. R. Castillo, Spectrochim. Acta, Part B, 2020, 169, 1–11.
3. I. V. Kubrakova, O. N. Grebneva-Balyuk, D. V. Pryazhnikov, M. S. Kiseleva and O. O. Efanova, J. Anal. Chem., 2023, 78, 1306–1319.
4. I. Kalomista, A. Keri and G. Galbacs, J. Anal. At. Spectrom., 2016, 31, 1112–1122.
5. A. R. Donovan, C. D. Adams, Y. Ma, C. Stephan, T. Eichholz and H. Shi, Chemosphere, 2018, 195, 531–541.
6. H. E. Pace, N. J. Rogers, C. Jarolimek, V. A. Coleman, C. P. Higgins and J. F. Ranville, Anal. Chem., 2011, 83, 9361–9369.
7. M. Cosmi, N. Gonzalez-Quiñonez, P. Tejerina Díaz, Á. Manteca, E. Blanco-González, J. Bettmer, M. Montes-Bayón and M. Corte-Rodriguez, J. Anal. At. Spectrom., 2021, 36, 2007–2016.
8. M. D. Montano, B. J. Majestic, A. K. Jamting, P. Westerhoff and J. F. Ranville, Anal. Chem., 2016, 88, 4733–4741.
9. S. A. Love, M. A. Maurer-Jones, J. W. Thompson, Y. S. Lin and C. L. Haynes, Annu. Rev. Anal. Chem., 2012, 5, 181–205.
10. T. Zhao, M. Fang, Z. Tang, X. Zhao and F. Wu, Sci. Total Environ., 2019, 692, 660–668.
11. M. Tharaud, A. P. Gondikas, M. F. Benedetti, F. von der Kammer, T. Hofmann and G. Cornelis, J. Anal. At. Spectrom., 2017, 32, 1400–1411.
12. S. Pradhan, J. Hedberg, E. Blomberg, S. Wold and I. Odnevall Wallinder, J. Nanoparticle Res., 2016, 18, 285.
13. S. Hellmann, T. Gil-Díaz, M. Böhm, D. Merten, S. Grangeon, F. Warmont, S. Unbehau, T. Sowoidnich and T. Schäfer, ACS Omega, 2024, 9(28), 30294–30307.
14. S. Dekkers, P. Krystek, R. J. B. Peters, D. P. K. Lankveld, B. G. H. Bokkers, P. H. Van Hoeven-Arentzen, H. Bouwmeester and A. G. Oomen, Nanotoxicology, 2011, 5, 393–405.
15. H. Nabeshi, T. Yoshikawa, K. Matsuyama, Y. Nakazato, K. Matsuo, A. Arimori, M. Isobe, S. Tochigi, S. Kondoh, T. Hirai, T. Akase, T. Yamashita, K. Yamashita, T. Yoshida, K. Nagano, Y. Abe, Y. Yoshioka, H. Kamada, T. Imazawa, N. Itoh, S. Nakagawa, T. Mayumi, S. ichi Tsunoda and Y. Tsutsumi, Biomaterials, 2011, 32, 2713–2724.
16. M. R. Go, S. H. Bae, H. J. Kim, J. Yu and S. J. Choi, Front. Microbiol., 2017, 8, 1–12.
17. E. Bolea-Fernandez, D. Leite, A. Rua-Ibarz, L. Balcaen, M. Aramendía, M. Resano and F. Vanhaecke, J. Anal. At. Spectrom., 2017, 32, 2140–2152.
18. J. Mosquera, I. García and L. M. Liz-Marzán, Acc. Chem. Res., 2018, 51, 2305–2313.
19. R. Thomas, in Practical Guide to ICP-MS, Marcel Dekker, New York, 1st edn, 2004, pp. 71–79..
20. M. Vázquez Peláez, J. M. Costa-Fernández and A. Sanz-Medel, J. Anal. At. Spectrom., 2002, 17, 950–957.
21. L. Hendriks, A. Gundlach-Graham, B. Hattendorf and D. Günther, J. Anal. At. Spectrom., 2017, 32, 548–561.
22. S. E. Szakas, R. Lancaster, R. Kaegi and A. Gundlach-Graham, Environ. Sci. Nano, 2022, 9, 1627–1638.
23. M. Resano, M. Aramendía, E. García-Ruiz, A. Bazo, E. Bolea-Fernandez and F. Vanhaecke, Chem. Sci., 2022, 13, 4436–4473.
24. I. Abad-Álvaro, E. Peña-Vázquez, E. Bolea, P. Bermejo-Barrera, J. R. Castillo and F. Laborda, Anal. Bioanal. Chem., 2016, 408, 5089–5097.
25. S. Theiner, K. Loehr, G. Koellensperger, L. Mueller and N. Jakubowski, J. Anal. At. Spectrom., 2020, 35, 1784–1813.
26. J. Liu, X. Wei, C. Wu, L. Zheng, M. Wang, M. Chen and J. Wang, Anal. Chim. Acta, 2023, 1254, 341114.
27. S. Fernández-Trujillo, M. Jiménez-Moreno, N. Rodríguez-Fariñas and R. C. Rodríguez Martín-Doimeadios, Anal. Bioanal. Chem., 2024, 416, 2657–2676.
28. F. Laborda, I. Abad-Álvaro, M. S. Jiménez and E. Bolea, Spectrochim. Acta Part B At. Spectrosc., 2023, 199, 106570.
29. S. Nunez, H. Goenaga, M. Yamanaka, T. Itagaki and S. Wilbur, Analysis of 10 nm gold nanoparticles using the high sensitivity of the Agilent 8900 ICP-QQQ, Agil. Appl. Note No. 5991-6944EN, 2016, pp. 1–6.
30. M. Yamanaka and S. Wilbur, Accurate Determination of TiO2 Nanoparticles in Complex Matrices using the Agilent 8900 ICP-QQQ, Agil. Appl. Note No. 5991-8358EN, 2017, pp. 1–7.
31. M. Yamanaka, T. Itagaki and S. Wilbur, High sensitivity analysis of SiO2 nanoparticles using the Agilent 8900 ICP-QQQ in MS/MS mode, Agil. Appl. Note No. 5991-6596EN, 2016, pp. 1–6.
32. T. E. Lockwood, R. Gonzalez De Vega and D. Clases, J. Anal. At. Spectrom., 2021, 36, 2536–2544.
33. G. Cornelis and M. Hassellöv, J. Anal. At. Spectrom., 2014, 29, 134–144.
34. G. Cornelis and S. Rauch, Anal. Bioanal. Chem., 2016, 408, 5075–5087.
35. G. Cornelis, Nanocount Software, https://blogg.slu.se/nanocount/, accessed 22 March 2024.
36. M. I. Chronakis, B. Meermann and M. von der Au, Anal. Bioanal. Chem., 2025, 417, 7–13.
37. F. Laborda, A. C. Gimenez-Ingalaturre, E. Bolea and J. R. Castillo, Spectrochim. Acta, Part B, 2019, 159, 105654.
38. L. A. Currie, Anal. Chem., 1968, 40, 586–593.
39. M. Corte Rodríguez, R. García Álvarez-Fernández, E. Blanco, J. Bettmer and M. Montes-Bayón, Anal. Chem., 2017, 89, 11491–11497.
40. P. García Cancela, N. González Quiñónez, M. Corte-Rodríguez, J. Bettmer, A. Manteca and M. Montes-Bayón, Metallomics, 2022, 14, mfac015.
41. M. D. Montaño, J. W. Olesik, A. G. Barber, K. Challis and J. F. Ranville, Anal. Bioanal. Chem., 2016, 408, 5053–5074.
42. F. Laborda, J. Jiménez-Lamana, E. Bolea and J. R. Castillo, J. Anal. At. Spectrom., 2013, 28, 1220–1232.
43. L. A. Currie, J. Radioanal. Nucl. Chem., 2008, 276, 285–297.
44. F. Laborda, E. Bolea and J. Jiménez-Lamana, Anal. Chem., 2014, 86, 2270–2278.
45. GeeksforGeeks, Python - Gaussian fit, https://www.geeksforgeeks.org/python/python-gaussian-fit/, accessed 22 March 2024.
46. Agilent Technologies, MassHunter Workstation Single Nanoparticle Application, Tokyo, 2019..
47. A. Laycock, N. J. Clark, R. Clough, R. Smith and R. D. Handy, Environ. Sci. Nano, 2022, 9, 420–453.
48. S. Hellmann, Application of single particle/single cell ICP-MS to geotechnical and biogeochemical suspensions, Friedrich Schiller University Jena, 2025,  DOI:10.17617/2.3654654..
49. S. Cuello-Nuñez, I. Abad-Álvaro, D. Bartczak, M. E. Del Castillo Busto, D. A. Ramsay, F. Pellegrino and H. Goenaga-Infante, J. Anal. At. Spectrom., 2020, 35, 1832–1839.

---