Simplistic Software for Analyzing Mass Spectra and Mixed Experimental-Theoretical Database for Identifying Poisonous and Explosive Substances

1. Introduction

The Chemical Weapons Convention [1], which entered into force in 1997, marked a breakthrough in a long-standing effort to end the production, storage, and eventual deployment of poisoning agents in a military setting. Despite its nearly universal adoption, multiple large-scale assaults involving chemical weapons have occurred in the decades after the adoption, most notably in Syria (before [2] and after [3] its accession to the convention) and Iraq [4]. In a concerning development, nerve combat agents, originally designed for indiscriminate large area use, have been employed in attempts on lives of individuals in urban environments. The most known case is the Tokyo subway sarin attack, performed in 1995 by the Aum Shinrikyo cult, that killed 13 and injured more than 6000 people [5,6]. In recent years, authoritarian regimes in Russia and North Korea [7,8] have made targeted attempts at using various poisons to assassinate dissidents and critics [9,10,11]. Thus, Russian democratic opposition leader Alexei Navalny [12,13] and a former Russian spy and double agent for British intelligence Sergei Skripal [14] were notoriously poisoned with the Novichok nerve agent, Ukrainian president Victor Yushchenko was poisoned during his presidential campaign of 2004 by the TCDD agent [15], and an exiled relative of North Korea’s supreme leader Kim Jong Un, Kim Jong–nam [16], was killed using VX nerve agent. Months after the attempt on Skripal, an unrelated British couple was poisoned with Novichok [17], apparently as a collateral of a Russian attack.

Although in the aforementioned high-profile cases the specific nerve agents were reliably identified, investigations into other apparent poisonings did not produce conclusive results on the nature of chemical agents used. In cases of Russian regime critics Pyotr Verzilov [18], Dmitry Bykov [19], Vladimir Kara–Murza Jr. [20], the latter being poisoned on two separate occasions, or in a recent chain of poisonings of dissident Russian journalists and activists after the outbreak of Russian aggression against Ukraine [21], the used substances were not definitively established, which might be due to delays in samples collection and their analysis.

A range of methods exists to identify the presence of chemical warfare agents in the laboratory or in the field. The most sensitive and informative of these are non-portable techniques: mass spectrometry (MS), nuclear magnetic resonance (NMR), and chromatographic methods, such as gas chromatography (GC) or high-pressure liquid chromatography (HPLC), coupled to MS [22,23,24,25,26]. In their review on the detection and destruction of chemical warfare agents, Kim et al. [27] provide numerous examples of MS techniques being used to identify organophosphorus nerve agents and other toxins at very low concentrations, in some cases in vivo. In most MS techniques, the molecules present in the sample undergo fragmentation upon ionization, which makes interpretation of mass spectra a cumbersome task even when dealing with a clean individual substance, increasing the likelihood of failure to identify a compound in the probe. In real-world forensic samples, often heavily contaminated and containing only traces of compounds, reliable identification becomes an exceedingly difficult task. Thus, a method to automatically identify poisons or other dangerous chemical compounds in mass spectra of impure samples is of great interest to a broad community of forensic experts, medical professionals, as well as independent sleuths. Since investigations are often conducted by individuals and teams with no technical education and at their own risk, we also note that a software piece to implement this method must be easy to install and operate without MS specialist knowledge.

Focusing on MS as the prime method to identify various species in experimental mixtures, we find ourselves with a wide selection of program tools for analyzing mass spectra. First of all, many producers of MS equipment provide accompanying software to be used with it. The MassHunter code by Agilent [28] is one such example. Secondly, the analysis software developed by the National Institute of Standards and Technology (NIST), such as the AMDIS (Automated Mass Spectrometry Deconvolution and Identification System) and MS Search [29,30,31,32,33] are commonly used. The drawback of these programs is that they are proprietary. As an alternative, there are also open-source software, such as the ProteoWizard [34], matchms [35,36], OpenMS/pyOpenMS [37,38], and FastEI [39]. However, most of these packages require both advanced user experience and proficiency in MS. Therefore, these software packages can be hard to use for non-experts.

Finding the reference spectra in the existing literature might also present a challenging problem. In the publicly accessible databases, such as those by The NIST Chemistry WebBook [40] or National Institute of Advanced Industrial Science and Technology (AIST) [41], experimental data for many substances are not present, for instance for the compounds described in the book by V.S. Mirzayanov [42]. There are some attempts to combine personal libraries of spectra, e.g, FederEI [43], a federated library matching framework for EI-MS. Another possible solution to this problem is to predict spectra from theory. Nowadays, various methods for such prediction exist. Among those are the machine-learning-based prediction algorithms, such as competitive fragmentation modeling (CFM) [44,45,46], rapid approximate subset-based spectra prediction (RASSP) [47], and neural electron-ionization mass spectrometry (NEIMS) [48]. In recent years, an algorithm to compute mass spectra by means of molecular dynamics (MD) simulations was proposed by S. Grimme [49]. This algorithm was used to predict the MS spectra, among others, of Tabun [50] and Novichok [51], experimental work therewith being greatly hindered by the inherent danger.

To address the outlined difficulties, we present a simple-to-use software package, ToxicMassSceptic, for the analysis of mass spectra, together with a database compiled from both MS experiments and theoretical computations, as well as the workflow for producing the theoretical mass spectra. We do not aim to outperform existing identification methods and libraries but rather to provide a simple and robust tool for preliminary substance identification that can aid low-budget analytical laboratories and civil investigators. The article has the following structure. First, in Section 2, we introduce the methodology: the structure and sources of the database, the digital formats of the data, and algorithms and workflows to compute and assign mass spectra, including the spectral similarity metrics. Secondly, we discuss the theoretical computation of mass spectra and demonstrate applications of the methodology in Section 3. Finally, conclusions are outlined in Section 4.

2. Methods

2.1. Mass-Spectroscopic Database

2.1.1. Database Structure and File Formats

Our database has to be easy to extend even by inexperienced users. Therefore, we store it as a set of nested directories with the structure shown in Figure 1. The top-level directory (“database”) contains the subdirectories that name the class of substances (“class #1”, “class #2”, etc.). Each of the subdirectories (“substance #1”, “substance #2”, etc.) contains folders with data on the specific substance. The recommended naming of these folders is “[Brutto chemical formula in the Hill notation]_[common name of the substance].” For every substance, the “ref.ms” file is required, which contains the reference mass spectrum of the given compound. It is optional but strongly suggested to supplement an entry with a file “INFO.txt” that contains information about the substance, e.g., common names, molar mass, links to substance Wikipedia and/or PubChem webpage, etc.

The classes of substances in the presented database and the number of entries in each are shown in Table 1. While the classification of substances is almost always self-explanatory, assuming their separation into different chemical weapon agent types (blister agents, blood agents, chocking agents, lachrymators, and nerve agents) [52], environmental pollutants (such as polycyclic aromatic hydrocarbons (PAHs) [53,54], per- and polyfluoroalkyl substances (PFAS), polychlorinated biphenyls (PCBs), and dioxines [55,56]), a separate category (miscellaneous) had to be made to store different substances that did not fit into this arguably rigid framework.

The reference spectra of the molecules in the database (files “ref.ms”) are formatted as two-column text files with pairs of numbers ($x,y$) in rows, where x is the integer mass-over-charge ($m/z$) position of the ion and y is the normalized intensity of the given ionic fragment in the MS; this format is usually denoted with an .xy file extension. The spectra in the “ref.ms” files have different normalization and are to be treated as not normalized, while normalization happens during the runtime. For a molecule with a spectrum of N fragment ions $\{(x_1,y_1),(x_2,y_2),\dots ,(x_N,y_N)\}$, the intensities are normalized such that

$${\displaystyle \sum _{i=1}^N}{y}_i=100\%.$$

(1)

2.1.2. Sources of Experimental Mass Spectra

Our database of molecular species borrowed mainly from the following sources: The NIST Chemistry WebBook [40], Spectral Database for Organic Compounds SDBS [41] organized by the AIST, Japan, and University of Rhode Island Explosives Database [57]. Since the Chemistry WebBook removed the option to download numerical MS data, most of the information from this database was extracted through manual digitizing the graphs (for details of this procedure, see ESI). The spectra for the two Novichok species, A-230 and A-232, were digitized from Ref. [58] using WebPlotDigitizer software [59].

2.1.3. Sources of Theoretical Mass Spectra

Theoretical mass spectra were computed using the workflow shown in Figure 2. All quantum chemical calculations, including conformational search and the MS calculation, were done with the GFN2-xTB method [60] as implemented in the xTB software [61], version 6.6.1. First, the initial molecular structure, obtained either from the NIST Chemistry WebBook, PubChem, or drawn in Jmol [62], was optimized with the xTB software. Then, conformational search was performed for this structure using CREST (version 2.12) [63,64], except for conformationally-rigid molecules. Subsequently, two augmented Born–Oppenheimer molecular dynamics (aBOMD) program packages were applied to calculate theoretical mass spectrum of the lowest energy conformer: QCxMS (version 5.2.1) [49,65,66], an original approach by S. Grimme, and DissMD, a software [67,68,69] based on the same idea. A detailed comparison of those approaches can be found in Section 3. Finally, the spectra obtained by the two theoretical approaches described above were combined as arithmetic means.

In QCxMS, the default settings were applied. The molecules were ionized by electron ionization (EI) with kinetic energy of electrons equal to 70 eV. The spectra were then collected by PlotMS (version 6.1.0). Since DissMD only simulates laser ionization, the ionization of molecules was modeled with an extreme ultraviolet (XUV) photon of 70 eV energy. In both QCxMS and DissMD calculations, the GFN2-xTB method was used to provide the potential energy surfaces for the aBOMD simulations, as this method was shown to be sufficiently accurate and computationally feasible for the mass spectra prediction [65,70].

2.2. Mass-Spectra Assigning Algorithm

2.2.1. Window-Function Based Assignment

The assignment was based on the assumption that there might be more than one species in the MS, which can be the case if the mixtures were not properly separated by chromatography or an alternative technique applied before the MS analysis. Therefore, the procedure involves finding only the relevant peaks in the tested spectrum to be compared with the reference database. For this, the window-based metrics were employed as described in more detail in the following.

Let us assume that we are interested in the possibility of species A with known reference spectrum of $N\left(A\right)$ peaks $\{(x_1,y_1),(x_2,y_2),\dots ,(x_{N\left(A\right)},y_{N\left(A\right)})\}$ to be present in the mixture. Intensities $y_i$ can be represented as an $N\left(A\right)$-dimensional vector $\mathbf{y}\left(A\right)=(y_1,y_2,\dots ,y_{N\left(A\right)})$. Note that we require all intensities to be positive ($y_i>0$ for $i=1,\dots ,N\left(A\right)$) and normalized to 100% as seen from Eq. (1). To make the comparison, we need to reduce the experimental dataset to an analogous $N\left(A\right)$-dimensional vector of experimental intensities $\tilde{\mathbf{y}}\left(B\right|A)=({\tilde{y}}_1,{\tilde{y}}_2,\dots {\tilde{y}}_{N\left(A\right)})$, where ${\tilde{y}}_i$ is the spectral intensity around $x_i=m_i/z_i$ in experimentally-measured MS $I\left(x\right)$ of unknown species or mixture B. To that end, we integrate the raw experimental MS $I\left(x\right)$ with a window function $w\left(x\right|x_i)$ for a given position $x_i=m_i/z_i$ and obtain non-normalized intensities $(Y_1,Y_2,\dots ,Y_{N\left(A\right)})$ as

$$Y_i={\displaystyle \underset{0}{\overset{+\infty }{\int }}}I\left(x\right)·w\left(x\right|x_i)dx,$$

(2)

where $w\left(x\right|x_i)$ is nonzero only in the vicinity of $x_i$. This mathematical operation essentially sums up the spectral intensity near an expected position $x_i$ into a single value. Applying this transformation to every peak i in the reference spectrum $I\left(x\right)$ and subsequently normalizing resulting values $Y_i$ such that

$${\tilde{y}}_i=100\%\times {\displaystyle \frac{Y_i}{{\displaystyle \sum _{j=1}^{N\left(A\right)}}{Y}_j}},$$

(3)

we obtain experimental intensities ${\tilde{y}}_i$ at the discretized positions $x_i=m_i/z_i$ of the reference dataset.

Alternatively, if the experimental MS is presented in the form of discrete peaks, the integration procedure is replaced by the summation, namely

$$Y_i={\displaystyle \sum _{k=1}^M}{I}_k·w\left(x_k\right|x_i)dx,$$

(4)

where index k runs over all M peaks with intensities $I\left(x_k\right)=I_k$ identified in the experimental MS by the spectrometer’s software.

In our program code, we implemented two types of window functions $w\left(x\right|x_i)$: A rectangular window,

$$w\left(x\right|x_i)=\left\{\begin{array}{cc}1\hfill & ,|x-x_i|\le \sigma /2,\hfill \\ 0\hfill & ,|x-x_i|>\sigma /2,\hfill \end{array}\right.$$

(5)

and Gaussian window

$$w\left(x\right|x_i)=exp\left(-{\displaystyle \frac{{(x-x_i)}^2}{2{\sigma }^2}}\right),$$

(6)

where $\sigma$ is the width of the given window in $m/z$ units. By default, the Gaussian window with $\sigma =1/2$ is employed.

2.2.2. Assignment Metric

After defining the window-based reduction scheme of experimental data, we can discuss the route to identifying chemical species in our spectrum. To that end, we rely on a metametric, which is composed of several deterministic metrics. Thus, the simplest metric $N\left(B\right|A)$ that can be defined for a given reference spectrum A is the number of lines present in both A and B. It reads

$$N\left(B\right|A)={\displaystyle \sum _{i=1}^{N\left(A\right)}}\theta ({\tilde{y}}_i-c),$$

(7)

where $c>0$ is a small threshold (in our case, $c={10}^{-15}$) for numerical comparison of real numbers and $\theta \left(x\right)$ is the Heaviside step function of the form

$$\theta \left(x\right)=\left\{\begin{array}{cc}1,\hfill & x>0,\hfill \\ 0,\hfill & x\le 0\hfill \end{array}\right..$$

(8)

The expression in Eq. (7) can be normalized by the total number of lines in the reference spectrum $N\left(A\right)$ to produce the relative number of lines, i.e.,

$$P\left(B\right|A)={\displaystyle \frac{N\left(B\right|A)}{N\left(A\right)}}.$$

(9)

More sophisticated metrics should also account for the distribution of fragment intensities. For this purpose, two sets of normalized values $\mathbf{y}\left(A\right)$ and $\tilde{\mathbf{y}}\left(B\right|A)$ can be treated as probability distributions. Thus, standard statistical distances for probability distributions can be employed. We chose four such measures: Kullback–Leibler divergence ($D_{\mathrm{KL}}$) [71], Bhattacharyya distance ($D_{\mathrm{B}}$) [72], Hellinger distance ($D_{\mathrm{H}}$) [73], and cosine distance ($D_{\mathrm{C}}$). In our case of two spectra, A and B, these four measures are given as [74,75]

$$D_{\mathrm{KL}}\left(B\right|A)={\displaystyle \sum _{i=1}^{N\left(A\right)}}{\tilde{y}}_i·ln\left({\displaystyle \frac{{\tilde{y}}_i}{y_i}}\right),$$

(10)

$$D_{\mathrm{B}}\left(B\right|A)=-ln\left(BC\left(B\right|A)\right),$$

(11)

$$D_{\mathrm{H}}\left(B\right|A)=\sqrt{1-BC\left(B\right|A)},$$

(12)

$$D_{\mathrm{C}}\left(B\right|A)=1-{\displaystyle \frac{{\displaystyle \sum _{i=1}^{N\left(A\right)}}{y}_i{\tilde{y}}_i}{\sqrt{\left({\displaystyle \sum _{i=1}^{N\left(A\right)}}{\tilde{y}}_i^2\right)·\left({\displaystyle \sum _{i=1}^{N\left(A\right)}}{y}_i^2\right)}}},$$

(13)

respectively. In Eqs. (11) and (12), $BC$ is the so-called Bhattacharyya dimensionless coefficient [72,76] given by

$$BC\left(B\right|A)={\displaystyle \frac{1}{100\%}}{\displaystyle \sum _{i=1}^{N\left(A\right)}}\sqrt{{\tilde{y}}_i·y_i}.$$

(14)

Here, the division by 100% is motivated by the fact of $BC$ being defined for probability distributions normalized to 1. The three chosen measures of similarities for probability distributions from Eqs. (9)–(12) require that components of the vector $\tilde{\mathbf{y}}\left(B\right|A)$ are non-negative. Note that Eqs. (10)–(12) are undefined for $N\left(B\right|A)=0$, which corresponds to the case of the species not being present in the spectrum.

The combined metametric is then constructed from Eqs. (9)–(13) such that

$$D_{\mathrm{meta}}\left(B\right|A)={\displaystyle \frac{1}{P\left(B\right|A)\times {\displaystyle \sum _{j=1}^{N\left(A\right)}}{Y}_j}}\times \left({\displaystyle \frac{D_{\mathrm{KL}}\left(B\right|A)}{{\varsigma }_{\mathrm{KL}}}}+{\displaystyle \frac{D_{\mathrm{B}}\left(B\right|A)}{{\varsigma }_{\mathrm{B}}}}+{\displaystyle \frac{D_{\mathrm{H}}\left(B\right|A)}{{\varsigma }_{\mathrm{H}}}}+{\displaystyle \frac{D_{\mathrm{C}}\left(B\right|A)}{{\varsigma }_{\mathrm{C}}}}\right),$$

(15)

where $Y_j$ is the non-normalized experimental intensity given by Eqs. (2) or (4) and ${\varsigma }_X$ is the standard deviation of the given metric $X=$ KL, B, H, and C, computed over the whole available dataset as

$${\varsigma }_X=\sqrt{\langle D_X^2\rangle -{\langle D_X\rangle }^2}=\sqrt{{\displaystyle \frac{1}{N_{\mathrm{d}}}}\sum _a{D}_X^2\left(B\right|A)-{\left({\displaystyle \frac{1}{N_{\mathrm{d}}}}{\sum }_a{D}_X^2\left(B\right|A)\right)}^2},$$

(16)

where index a runs over all spectra in the database and $N_{\mathrm{d}}$ is the number of such spectra. The value of $D_{\mathrm{meta}}\left(B\right|A)$ from Eq. (15) tends to zero if the two spectra A and B are similar and increases with the growing dissimilarity of the experimental spectrum from the reference. Although Bhattacharyya and Hellinger distances provide the same relative ranking of substances, it can be advantageous to use both in the metametric, as they might have different sensitivity at different values of the Bhattacharyya dimensionless coefficient $BC$.

2.2.3. Background Removal Algorithm

Experimentally measured spectra can contain signals from the background. This may result in empty areas of a spectrum producing negative intensities when using Eqs. (2) and (4). To avoid that, basic filtering of the experimental MS signal $I\left(x\right)$ can be performed. The simplest and most robust approach is probably a visual determination of the noise threshold level $I_{\mathrm{thr}}$, and setting all the values $I\left(x\right)\le I_{\mathrm{thr}}$ to zero. However, a crude automatic routine can also be designed (for example, see Ref. [77]) assuming that non-zero peaks occupy only a minor part of the spectrum in all available $m/z$ ranges and that the baseline signal is $I=0$. To that end, we represent a spectrum in a discretized form with lines $I_1,I_2,\dots ,I_M$. Then, the following procedure can be employed.

• Calculate the standard deviation of $I\left(x\right)$ from baseline ($I=0$) as ${\mathrm{SD}}_0=\sqrt{{\displaystyle \frac{1}{M}}{\sum }_{k=1}^M{I}_k^2}$.
• Consider only values $I_k<q·{\mathrm{SD}}_0$, with $q\ge 1$ being an arbitrary selectivity coefficient, forming a new set $I_1^{\left(1\right)},I_2^{\left(1\right)},\dots ,I_{M_1}^{\left(1\right)}$, where the upper index “(1)” indicates the iteration number and $M_1\le M$ is the number of elements in the new set.
• Calculate the new standard deviation as ${\mathrm{SD}}_1=\sqrt{{\displaystyle \frac{1}{M_1}}{\sum }_{k=1}^{M_1}{\left(I_k^{\left(1\right)}\right)}^2}$.
• Repeat steps 2 and 3 until the number of elements in the set remains constant or a maximum number of iterations p is reached.
• Set values of the original mass spectrum below the final threshold $q·{\mathrm{SD}}_p$ to zero.

This automatic background removal procedure is implemented in our program code, with the default number of steps $p=3$ and selectivity coefficient $q=1.5$.

2.3. Software

The program code called ToxicMassSceptic is written in Python version 3.8 for the Linux, MacOS, and MS Windows operational systems, distributed under an open source Apache License version 2.0 [78], and is managed using the version control system GIT [79] by the provider GitLab [80]. The source code is available under Ref. [81]. The list of program requirements includes Python packages such as Numpy [82] and Matplotlib [83]. The code has a clear version number and is accompanied by two types of documentation: (i) a README file in the Markdown format outlining external dependencies, package structure as well as the installation procedure and (ii) an automatically generated Doxygen [84] code documentation describing all constituting objects and functions. The package-management system PIP3 [85] governs the installation procedure. The code is aimed to be fully unit-tested. To that end, the package Unittest [86] is employed. The current code design enables the use of our program as an external Python library as well as through a command line interface.

The flowchart of the ToxicMassSceptic work and usage is given in Figure 3. First, the user needs to provide a spectrum, which can then be passed, by request, through the background removal procedure described in Section 2.2.3. Then, the database is loaded, and the comparison of the unknown spectrum with the database entry begins. During this step, the four metrics described in Section 2.2.2 are computed for each substance. After all the metrics are known, the metametric from Eq. (15) for each database substance is computed, as it requires a spread of each metric throughout every database entry as seen from Eq. (16). Finally, the database entries are sorted by the metametric value, and the best matching substance is given.

2.4. Statistical Analysis of Results

Let us assume that the user is interested in testing $N_{\mathrm{trials}}$ number of different mixtures B. Each such ith mixture $B_i$ contains a compound $A_{\mathrm{true}}^i$, which is also present in the database. Furthermore, we assume that for each sample $B_i$, the top-K matching candidates ${\mathbf{A}}^i=\{A_1^i,A_2^i,\dots ,A_K^i\}$ are suggested by our algorithm based on the metrics introduced above in Section 2.2.2. Here, each set ${\mathbf{A}}^i$ is sorted in descending order such that its first element is the most probable match. Therefore, index j denotes the rank of compound $A_j^i$, i.e., $j=R\left(A_j^i\right)$, with lower ranks being preferable. Then, the following scores can be introduced to assess the performance of our algorithm.

• Top-K accuracy (also known as Hit rate at rank K), which is equal to the number of trials with the correctly identified compound being present in top-best K candidates $N_{\mathrm{in}\mathrm{top}-K}$ divided by the total number of trials $N_{\mathrm{trials}}$ and multiplied by 100%, i.e.,

  $$\mathrm{top}-K\mathrm{accuracy}={\displaystyle \frac{N_{\mathrm{in}\mathrm{top}-K}}{N_{\mathrm{trials}}}}\times 100\%.$$

  (17)
• Mean reciprocal rank (MRR), defined as

  $$\mathrm{MRR}={\displaystyle \frac{1}{N_{\mathrm{trials}}}}{\displaystyle \sum _{i=1}^{N_{\mathrm{trials}}}}{\displaystyle \frac{1}{R\left(A_{\mathrm{true}}^i\right)}}\times 100\%,$$

  (18)

  where $R\left(A_{\mathrm{true}}^i\right)$ is the rank of the correctly identified compound $A_{\mathrm{true}}^i$ in trial i.
• Mean rank (MR), defined as

  $$\mathrm{MR}={\displaystyle \frac{1}{N_{\mathrm{trials}}}}{\displaystyle \sum _{i=1}^{N_{\mathrm{trials}}}}R\left(A_{\mathrm{true}}^i\right).$$

  (19)

The top-K score from Eq. (17) shows how often the correctly identified compound was present in the K most probable candidates predicted by the program code, whereas MRR from Eq. (18) evaluates the ability of the code to assign low ranks to relevant chemical compounds. In the case of an ideal assignment, when correct compounds always occupy the very top of the suggestion list, both scores are equal to 100%. The MR score from Eq. (19) is closely related to MRR, but is equal to or greater than 1.0 and tends to 1.0 for better performing recommendation systems.

3. Results and Discussion

3.1. Mass-Spectra Prediction Workflow

Predicted mass spectra presented in this work were computed using either QCxMS or DissMD. The latter is a part of the PyRAMD package [67,87,88]. Both algorithms employ Born–Oppenheimer molecular dynamics (BOMD), as proposed by S. Grimme in his seminal paper [49]. Before discussing our results, we first compare those two approaches.

A graphical representation of an aBOMD-based theoretical workflow for an MS spectrum prediction is depicted in Figure 4. First, multiple molecular geometries are generated, representing the gaseous ensemble of molecules in the spectrometer. Those structures are then used as initial points to start BOMD dynamics for ions. To include the electronic excitation effects, the BOMD dynamics are perturbed (or augmented) by the kinetic energy influx from an external energy reservoir, producing aBOMD trajectory. This energy, referred to as the internal excess energy (IEE), and the ion charge are ascribed according to the ionization procedure. If, upon the aBOMD trajectory propagation, a dissociation of the molecule is detected, the parent ion trajectory is stopped, and new aBOMD trajectories for the products are initiated by sharing the charge and IEE of the parent ion between fragments. Then, these trajectories of the daughter ions are propagated further. Finally, the mass spectra are computed from the ensemble of MD trajectories by counting the final products.

Despite this scheme’s general simplicity, a few crucial components in the algorithm define the simulation behavior. The QCxMS and DissMD use two completely different approaches to generate initial conditions. In the QCxMS, the thermostated MD of the neutral molecule is performed to sample the initial structures and their velocities. In the DissMD, the simplified Wigner sampling [88,89] approach from a user-provided geometry is used, which, in principle, can include some of the nuclear quantum effects [90] for the lighter nuclei such as hydrogens. Furthermore, these two approaches also differ greatly in the ionization procedure and the assignment of the IEE. In QCxMS, an arbitrary Poisson-like distribution is employed [65,91]

$$P\left(\mathrm{IEE}\right)={\displaystyle \frac{exp\left(c·\mathrm{IEE}·(1+ln\left(b/(c·\mathrm{IEE})\right)-b\right)}{a·\mathrm{IEE}+1}},$$

(20)

where $P\left(\mathrm{IEE}\right)$ is the probability of the ion to have the value of $\mathrm{IEE}$ upon ionization, whereas $a=0.2$ eV, $b=1$ eV, and $c=1/N_{\mathrm{ve}}$ are pre-defined parameters with $N_{\mathrm{ve}}$ being the number of valence electrons in the system. In the DissMD, however, an approach based on the electronic density of the states is used. Upon applying the maximum entropy principle and energy conservation to molecular ionization, one arrives at the following distribution [68]:

$$P\left(\mathrm{IEE}\right)=\mathrm{DoS}\left(\mathrm{IEE}\right)·{\left(E_{\mathrm{i}}-\mathrm{IP}-\mathrm{IEE}\right)}^{{\displaystyle \frac{N_{\mathrm{f}}}{2}}-1},$$

(21)

where $\mathrm{DoS}\left(\mathrm{IEE}\right)$ is the electronic density of states of the ion, $E_{\mathrm{i}}$ is the total energy of the ionization event, $\mathrm{IP}$ is the sum of ionization potentials to reach a given ionization state, and $N_{\mathrm{f}}$ is the number of degrees of freedom for the leaving particles. For the photoionization, which is the only available case in DissMD, $E_{\mathrm{i}}=mh\nu$ and $N_{\mathrm{f}}=3·N_{\mathrm{re}}$. In these expressions, m is the number of absorbed photons, $h=6.626\times {10}^{-34}$ J$·$s is the Planck constant, $\nu$ is the photons’ frequency, and $N_{\mathrm{re}}$ is the number of electrons removed upon ionization ($N_{\mathrm{f}}=3$ for single ionization, $N_{\mathrm{f}}=6$ for double ionization, etc.). Note, however, that Eq. (21) can still be applied for the electron impact ionization. In this case, $E_{\mathrm{i}}$ is the kinetic energy of the electrons and $N_{\mathrm{f}}$ is set to $3·(N_{\mathrm{re}}+1)$ to account for the leaving ionizing particle’s degree of freedom. Unlike in the first version of the software, in which the explicitly computed excited states were used to obtain the electronic density of states [87], the current version of the DissMD uses a simplified heuristic model based on the Van-der-Waals volume and surface to approximate $\mathrm{DoS}\left(\mathrm{IEE}\right)\propto {\mathrm{IEE}}^n$ as a power function with a single parameter n. In this case, Eq. (21) reduces to a beta-distribution [68].

The third crucial component of the simulation is the rate of internal conversion (IC), showing how fast the $\mathrm{IEE}$ decays into nuclear motions. For this purpose, the QCxMS uses the energy-gap law in the form [65]

$$k_{\mathrm{IC}}^{-1}={\displaystyle \sum _{j>i}^M}{\displaystyle \frac{k_{\mathrm{h}}}{N_{\mathrm{ve}}}}exp\left(\alpha ({\epsilon }_i-{\epsilon }_j)\right),$$

(22)

where $k_{\mathrm{h}}=2$ ps and $\alpha =0.5$ eV⁻¹ are constrants, ${\epsilon }_i$ is the energy of an i-th orbital, and M is the total number of orbitals. Contrary to that, in the DissMD a classical model of hot electrons with kinetic energy of $\mathrm{IEE}$ colliding with motionless nuclei is employed. In the DissMD prototype, a similar algorithm, based on an idea of electron–nuclear collision-induced IC, was used to compute the IC rates using the atomic electronic densities through the plasma frequency estimated from atomic charges [68]. However, in the newer code, it was replaced with a simplified model for the rate of such collisions is given as [69]

$$k_{\mathrm{IC}}=\kappa {\displaystyle \frac{\sqrt{m_{\mathrm{e}}\mathrm{IEE}}}{m_{\mathrm{amu}}(L_0+L_{\mathrm{mol}})}}{N}_{\mathrm{n}}{N}_{\mathrm{e}},$$

(23)

where $N_{\mathrm{e}}$ and $N_{\mathrm{n}}$ are the total number of electrons and nuclei in the ion, respectively, $m_{\mathrm{e}}$ is the electron mass, $m_{\mathrm{amu}}$ is the atomic mass unit (dalton), $L_{\mathrm{mol}}$ is the molecular length (atomic-charge-product-weighted sum of all chemical bonds, determined from the covalent radii of atoms), $L_0=5$ Å is the regularizing parameter, and $\kappa \approx 1.28$ is the fitted parameter based on the available experimental data [69].

When the dissociation is detected, the QCxMS and DissMD again proceed in a different fashion. The DissMD follows a direct route: Upon the detection of dissociation of ion ${\mathrm{M}}^{q+}$ into fragments $\mathrm{A}$ and $\mathrm{B}$ it calculates the energies of several channels

$${\mathrm{M}}^{q+}\to {\mathrm{A}}^{q_{\mathrm{A}}+}+{\mathrm{B}}^{q_{\mathrm{B}}+},$$

(24)

that satisfy the charge conservation $q_{\mathrm{A}}+q_{\mathrm{B}}=q$. Upon dissociation, the channels with non-negative kinetic energy release (KER) are assigned a probability proportional to this KER value. Subsequently, one of these channels is chosen according to those probabilities. This leads to a speedup in the calculation, as the neutral fragments are not propagated. However, this approach requires a larger number of trajectories to be computed. In the QCxMS, a concept of statistical charge, or statistical weighing, is used. In this approach, the MD is carried out for all fragments, but their associated intensities depend on the weight, which is determined as [65]

$$C_i={\displaystyle \frac{exp\left(-\frac{{\mathrm{IP}}_j}{k_{\mathrm{B}}T}\right)}{{\sum }_{j}exp\left(-\frac{{\mathrm{IP}}_j}{k_{\mathrm{B}}T}\right)}}$$

(25)

with indices i and j running over the number of fragments, ${\mathrm{IP}}_j$ being the ionization potential of a given fragment, $k_{\mathrm{B}}$ being Boltzmann constant, and $T=\mathrm{KE}/\left(3k_{\mathrm{B}}{N}_{\mathrm{n}}\right)$ being the instant temperature of nuclei, as computed from their kinetic energy (KE). With these fragment weights, it is also possible to directly apply the isotopic distribution in the post-analysis, while the DissMD requires running simulations with different isotopes.

To demonstrate the predictive capabilities of the aBOMD-based approach for computing the mass spectra, we took four molecules, for which we had the available spectra: methanol (${\mathrm{CH}}_3\mathrm{OH}$), novichok A-230 (${\mathrm{C}}_7{\mathrm{H}}_{16}{\mathrm{FN}}_2\mathrm{OP}$), o-chlorophenoxyacetic acid (${\mathrm{C}}_8{\mathrm{H}}_7{\mathrm{ClO}}_3$), and vinclozolin (${\mathrm{C}}_{12}{\mathrm{H}}_9{\mathrm{Cl}}_2{\mathrm{NO}}_3$). Structures of the most stable conformers of these molecules, according to CREST, can be found in Figure 5. As a metric to judge the similarity between spectra, we chose the number of peaks from the reference spectrum from Eq. (7), the Kullback–Leibler divergence given in Eq. (10), and the Bhattacharyya distance from Eq. (11).

The results of our comparison are demonstrated in Figure 6, Figure 7, and in Table 2. It is clear that the QCxMS, as the software specifically designed for EIMS predictions, outperforms DissMD. Nevertheless, in three out of four cases, DissMD provided extra fragments, which were missing in the QCxMS predictions. In all cases, the combination of the both methods allowed us to cover more than 80% of lines from experimental spectra. However, the relative intensities of the peaks are not always perfect, which can be a result of wrong ionization conditions in the simulations. Nevertheless, we can confirm the conclusions from previous studies in Refs. [50,51], stating that it is possible to use theoretically predicted mass spectra for the assignment of species with absent experimental reference spectra.

However, we would also claim that new software is probably due to development that would take the best algorithmic solutions from the QCxMS and DissMD. For the ionization stage, it makes more sense to assign the IEE from a physically sound model from Eq. (21). For computing the IC rate, one might use a better model of the electron-phonon coupling. One such possibility is demonstrated in Refs. [92,93], where the rate is calculated based on the Fermi–Dirac distribution and orbital overlaps for the two consecutive MD steps. For the treatment of dissociation, the QCxMS approach appears more suitable. However, instead of using the heuristically defined weights from Eq. (25), it would make more sense to use a modified version of the model introduced in Ref. [94], as it takes into account not only the ionization energies of fragments, but also the electron affinities, and the dissociation energies.

3.2. Performance Tests with Simulated Data

The ToxicMassSceptic features are subject to unit tests, ensuring the code works as expected. One of the production tests trials the performance of the code in the presence of noise and additional substances. Here, we perform the testing based on the undecayed substances in our database, simulating mixtures taken directly from the environment, rather than from biological samples. This is due to the fact that the analysis of biological substances usually requires the use of liquid chromatography and searching for metabolites, which can be known only from in vitro studies (see e.g., Ref. [25]). The biochemical degradation pathways of such compounds are highly unlikely to be found in publicly available literature in sufficient amounts to train empirical models. The theoretical prediction of those products from first principles is doubtful due to the sheer complexity of the problem. Therefore, without a proper database, we do not anticipate using ToxicMassSceptic directly for biological samples.

We model the species’ spectra with Gaussian-shaped peaks with randomly chosen standard deviation, that is, in the range between 0.05 and 0.1 $m/z$. We take the mass spectrum of a randomly chosen species from the database and generate a spectrum in the $m/z$ range from 0 to 500 with 2000 points. Then, we add a background that consists of two components. First, the signal of the substance is mixed with a spectrum composed of signals from benzene (${\mathrm{C}}_6{\mathrm{H}}_6$), oxygen (${\mathrm{O}}_2$), nitrogen (${\mathrm{N}}_2$), carbon dioxide (${\mathrm{CO}}_2$), and farnesene (${\mathrm{C}}_{15}{\mathrm{H}}_{24}$), one of sesquiterpenes. The relative amounts of the background species are randomly chosen between 0.1 and 0.2. Then, a random uniformly distributed noise is added on top of that with a signal-to-noise (S/N) level randomly chosen from the interval between S/N=100 and S/N=1000. Then, this generated spectrum is passed through our assignment algorithm, including the background removal and the rating of the actual compound, which is stored. The mean rating of the spectra upon multiple trials should not exceed an MR (Eq. 19) threshold, which, in our case, is set to 5. The current version of the software routinely passes this test.

To further demonstrate the performance of our code and compare different metrics, we carried out assignments of 500 randomly generated spectra. To that end, we modified the settings described above by lowering the allowed signal-to-noise level to $5\le \mathrm{S}/\mathrm{N}\le 100$, and additionally allowing peak intensities to vary by $\pm 50$% and their positions to be shifted by $\pm 0.2$ $m/z$. The assignment was repeated 48 times, leading to 24,000 trials in total and allowing us to compute the mean values and standard deviations for statistical parameters from Eqs. (17)–(19). The results of this analysis are shown in Table 3. As can be seen, the worse top-1 result is obtained using the cosine distance $D_{\mathrm{C}}$, reaching an accuracy level of only about 30%. The performance of other metrics is much higher and varies from about 55 to 91%. Similar trends are observed for the MRR and MR scores. The use of the proposed metametric $D_{\mathrm{meta}}$ was found to produce results of the highest quality in all cases.

3.3. Performance Test with Experimental Noisy Dataset

As an example of the mass spectra with noisy background, we took the strong-field-induced mass spectra of a tree-ring PAH fluorene (${\mathrm{C}}_{13}{\mathrm{H}}_{10}$), which are openly available from Ref. [95]. Since fluorene is in the database, and the laser-induced fragmentation patterns look similar to those obtained with EI, we simply tested the identification of the species with the mass spectra obtained using different laser peak powers (from $1.5\times {10}^{13}$ to $6.8\times {10}^{13}$ W/cm²). In all of the cases, the automatic background removal was applied.

The background removal results are shown in Figure 8. As one can see, the background is indeed removed quite efficiently, leaving only the signals from the ion fragments. The cleaning in the range of higher masses is somewhat less effective, which is due to the overall background level increase, as clearly seen in a logarithmic plot. Nevertheless, such background removal was sufficient to identify fluorene in the case of all experimental spectra considered in this work. The results for the highest peak power spectrum are shown in Figure 9.

3.4. Performance Tests with an Experimental Dataset of Cleaned Spectra

The mass spectra of 64 substances were recorded using GC (HP6890, Agilent Technologies) coupled to a single quadrupolar MS (HP5972A or HP5973, Agilent Technologies) or with GC (Trace 1310) coupled to MS (TSQ Duo Triple Quadrupole, Thermo Scientific). Helium was used as a carrier gas, and the spectra were measured in the range of 50–500 $m/z$. The EI was used to ionize species with an electron KE of 70 eV. More details on the measurement parameters are available in ESI.

The experimental dataset consists of several classes of substances: Acid contaminants, chlorophenols, dioxins, PAHs, pesticides, and herbicides. For each of the compounds from this dataset, the reference spectrum was added to the database, and then ToxicMassSceptic was tested to provide the assignment results. We ranked the performance in each dataset using six scores: Top-1, top-3, top-5, and top-10 accuracies from Eq. (17), MRR from Eq. (18), and MR from Eq. (19). The results of the test are given in Table 4. As one can see, most of the species were correctly identified in the top-3 best-matched substances, and the correct compound was the best-matched one 60% of the time, on average. With that, we conclude that the current performance allows the identification of species in unknown samples.

3.5. Testing Theoretical Reference Against Cleaned Experimental Data

In the dataset used in Section 3.4, there were three dioxines: 1,2-Dichlorodibenzo-p-dioxin, 1,3-Dichlorodibenzo-p-dioxin, and 1,4-Dichlorodibenzo-p-dioxin. These compounds are suitable for testing the assignment of experimental spectra against theoretically predicted mass spectra. For that reason, we computed the theoretical mass spectra of these three structural isomers using the workflow shown in Figure 2. In addition to that, in Section 3.1, we calculated theoretical mass spectra for o-chlorophenoxyacetic acid and vinclozolin, which were also present in the same database.

Thus, we took these five substances to test their identification with the ToxicMassSceptic software. The resulting ranking of these theoretical spectra (R) against their experimental counterparts is given in Table 5 in columns $\mathrm{Threshold}=0\%$. As one can see, the results are acceptable. However, upon examination of the theoretical spectra, one can see that the number of reference lines ($N_{\mathrm{lines}}$) is much larger than usually available for experimental spectra taken from various databases (which is typically of the order of a few tens of data points). Therefore, we have tried to remove some of the fragments with lower intensities from the theoretical spectra to see the effect on the identification of substances. In particular, we removed every lower-intensity peak by setting a relative threshold with respect to the most intensive one. We tried two settings: thresholds of 1% and 5%, which drastically reduced the number of lines and had an effect on the prediction performance (see Table 5). With 1% threshold, the MR value for this set of five spectra was slightly lower than the 0% and 5% settings, which indicates that there is an optimal amount of lines to represent a species in the database, as too many or too little may lead to misidentification of the species. Therefore, we recommend removing the weak intensity fragments when using ToxicMassSceptic for predicting theoretical mass spectra, as this improves the identification probability $P\left(B\right|A)$ (Eq. 9). The importance of the latter can be seen from the definition of the metametric from Eq. (15).

4. Conclusions

In this paper, we have presented an algorithm and a computer program for identifying toxic and combat compounds using mass spectrometry, ToxicMassSceptic, that is easy to operate for nonprofessionals. An essential part of it is the database of substances, assembled from multiple different sources, most prominently from databases like the NIST Chemistry WebBook and the SDBS of AIST, as well as from quantum chemical modeling. The use of theoretically predicted mass spectra allowed us to obtain reference data for poisonous substances for which no publicly accessible data exist. According to our tests against simulated and experimental datasets, ToxicMassSceptic with the database can facilitate preliminary identification of possible traces of poisonous and explosive substances. However, the current approach implies that the best matching result is always given. This is due to the open problem of finding thresholds for the current definition of the metametric. Therefore, the identification results are always biased toward the available database. The preliminary analysis results should always invoke a manual inspection for a few of the best-matching substances, to check that the identified peaks are indeed present in the spectra. The final conclusions regarding substance identification should always be based on expert opinion and validated with other experimental methods, such as NMR or rotational spectroscopy [96].