Comparison of in-silico and experimental spectra of example molecules
Following the general workflow, we first tested the QCEIMS software on two trajectories for a simple molecule, 3-cyclobutene-1,2-dione (Figure 2). The observed fragment ions yielded an excellent weighted dot-product similarity score of 972 and a cosine similarity of 839. When analyzing the trajectories to show the fragmentation pathways, we found clear evidence of the mechanisms by which the three main product ions observed in the experimental mass spectrum were produced (m/z 82, 54, 26), i.e., molecular ion, a neutral loss of carbon monoxide [M-CO]+ and loss of another carbon monoxide to yield [M-2CO]+ (Figure 2). Trajectory 2 lasted only 402 fs until the maximum of three fragments per trajectory was achieved (set in the QCEIMS source code), while trajectory 1 lasted 656 fs, because the initial two fragments reached a stable state and did not fragment further for a long time. The QCEIMS predictions also agreed with mechanisms predicted by the heuristic rule-based commercial MassFrontier software, showing first an α-cleavage followed by a CO molecule loss. This simple example shows that QCEIMS can generate correct molecular fragments and predict reasonable reaction mechanisms.
Here we show six molecules (Figure 3a-3f) as examples for QCEIMS predicted spectra versus experimental library spectra (Table 1). These examples demonstrate that QCEIMS yields different prediction accuracies. The examples also show different degrees of molecular flexibility. For each molecule, spectra showed specific characteristics that are here explained in brief.
Table 1. Mass spectral similarities of QCEIMS simulations against experimental spectra for select compounds
Name
|
InChIKey (short)*
|
M.W.**
|
RBN
|
PHI
|
Dot
|
Cos
|
2,4-Dimethyl-oxetane
|
KPPWZEMUMPFHEX
|
86.07
|
2
|
2.64
|
414
|
729
|
2-Nonene
|
IICQZTQZQSBHBY
|
126.27
|
5
|
7.52
|
789
|
762
|
2-Propynyloxy Benzene
|
AIQRJSXKXVZCJO
|
132.06
|
0
|
1.17
|
379
|
426
|
Furan
|
YLQBMQCUIZJEEH
|
68.08
|
0
|
0.55
|
988
|
806
|
1,8-Nonadiene
|
VJHGSLHHMIELQD
|
124.25
|
6
|
7.05
|
163
|
713
|
Adamantane
|
ORILYTVJVMAKLC
|
136.13
|
0
|
1.18
|
883
|
678
|
* first 14-characters of full InChIKey; **M.W. is the molecular weight in Daltons (Da); RBN (rotatable bond number) and PHI (Kier flexibility index) are rigidity descriptors and Dot and Cos are mass spectral similarity scores.
2,4-dimethyl-oxetane (Figure 3a): With a weighted dot-product score of 417, this spectrum represents a low quality in-silico prediction. We need to clarify that, for simplicity, we only calculated the spectrum of cis-2,4-dimethyl-oxetane, while its reference spectrum in NIST 17 mass spectral library contains no stereochemistry information because neither EI-MS nor chromatography technology can easily differentiate diastereomers. The experimental spectrum showed a low-intensity [M]˙+ at m/z 86 and initial neutral losses of a methyl-group and water (m/z 71 and m/z 68). QCEIMS did not predict these initial losses. Indeed, the high number of experimental fragment ions suggest that this molecule splits readily along multiple reaction pathways, most likely through breaking the molecular ether-bonds that subsequently break into smaller fragments. The main fragment ions at m/z 42 and m/z 44 were correctly predicted by QCEIMS as C3H6+ and C2H4O+ but not by the rule-based software MassFrontier. This case suggests that quantum mechanics-based simulations can produce novel reaction pathways that are absent from rule-base software predictions.
2-Nonene (Figure 3b): The in-silico spectrum of 2-nonene was highly similar to the experimental spectrum with dot-product match of 789. The main fragment ion at m/z 55 and the [M]˙+ at m/z 126 were very well reproduced. However, ion abundances of [M-1]+, [M-2] +, [M-3] + and [M-4] + were overestimated. In QCEIMS, these ions resulted from loss of several atomic or molecular hydrogens, suggesting that these bonds were fragmented more easily under semiempirical methods [23] than under experimental conditions.
Aromatic systems (Figure 3c and 3d): Both 2-propynyloxy benzene and furan were aromatic oxygen-containing molecules with low PHI values (1.71 and 0.55, respectively). Although the presence of most fragment ions was correctly predicted by QCEIMS for both molecules, dot-product similarity scores were radically different with a dot-product of 379 for 2-propynyloxy benzene and a dot-product similarity of 988 for furan). For 2-propynyloxy benzene, this low matching score was caused by the absence of an experimental [M]˙+ at m/z 132 that was largely overestimated in the in-silico spectrum. The fragmentation base ion (at 100% intensity) at m/z 93 represents the stable phenol ion and a neutral loss of C3H3, while the experimentally observed fragment at m/z 95 was missed in the QCEIMS prediction. At the same time, the presence of the C3H3+ product ion at m/z 39 (and a neutral loss of a phenol moiety) was overestimated by the QCEIMS method. This result suggests that the QCEIMS method needs further optimization in predicting the correct assignment of cation stability and assignment of the molecule with the lowest ionization energy in the fragmentation process (Stevenson’s rule [28]).
1,8-nonadiene (Figure 3e): For this molecule, a great disagreement between the cosine similarity score of 713 and the weighted dot-product of 163 was observed. The weighted dot-product emphasizes high m/z ions that are penalized if missing in spectral matching. Again, QCEIMS overestimated the abundance of the molecular ion [M]˙+ and of several atomic or molecular hydrogens from it. In addition, QCEIMS underestimated a neutral methyl loss (to m/z 109) and a neutral loss of ethylene (to m/z 96). To capture all potential fragmentations in QCEIMS such as the missed ethylene loss, more accurate PES estimates are needed.
Adamantane (Figure 3f): Adamantane is a well-known inflexible molecule. Our QCIEMS simulations correctly predicted the structure of the m/z 79 product ion as protonated benzene, proved by an independent publication of an infrared multiphoton dissociation spectrum [29] and DFT computations [30]. In comparison, the rule-based MassFrontier generated less reasonable fragment molecules that included cyclopropyl-moieties. The QCEIMS results showed that the m/z 93 product ion is likely associated with both ortho- and para-protonated toluene, in accordance with infrared multiphoton dissociation spectrum results [29]. These instances highlight the ability of QCEIMS to predict non-obvious mechanisms, such as rearrangements from sp3 hybrid carbons to aromatic system.
Probing the QCEIMS parameter space
A number of parameters can be chosen in the QCEIMS software, including the number of trajectories, impact excess energy per atom and initial temperatures. Other parameters such as the type of energy distribution and maximum MD time were excluded because they were already optimized during the development of QCEIMS [8]. We used OM2 because other semiempirical methods had been shown previously to perform worse [8]. For each molecule we chose one conformer and performed QCEIMS simulations with different parameter settings. By repeating QCEIMS simulations 50 times, we confirmed that identical mass spectra were obtained when using the same conformer under the same parameter settings. We changed parameter settings for 2,4-dimethyl-oxetane, 2-nonene and adamantane.
(1) Number of trajectories (ntraj)
In molecular dynamics, different reaction trajectories must be explored to cover possible routes of independent fragmentations across the energy surface. Each trajectory requires computational time, and therefore, the number of trajectories should be as low as possible. However, it is not clear a priori how many trajectories sufficiently cover the chemical reaction space and allow convergence to a consensus spectrum. By default, the QCEIMS program automatically calculates the number of trajectories by multiplying the number of atoms by 25. We explored this default value ranging from 8 to 1000 trajectories per atom for the different molecules, yielding up to 15,000 trajectories in total (Figure 4a). For each of the three molecules, the difference between the best and the worst similarity score differed only by 10% or less. None of the three molecules had improved similarity scores with higher number of trajectories. Indeed, it appeared that increasing the number of trajectories might lead to slightly lower dot-product similarity scores as observed for 2-nonene and adamantane, possibly due to a higher contribution of rare fragmentation reactions that lead to low abundant fragment ions that negatively impact similarity to experimental spectra. We concluded that the default value of 25 trajectories per atom number in a molecule was reasonable.
(2) Impact excess energy per atom (ieeatm)
Next we tested the impact excess energy (IEE) that is introduced by the colliding electron in electron ionization as vibrational energy into the molecules. The default value (ieeatm) in QCEIMS software is set at 0.6 eV per atom on the basis of previous OM2 tests [31]. At the beginning of each molecular dynamics simulation the molecule is heated by increasing the atom velocities until the impact excess energy is converted to kinetic energy that leads to bond fragmentation. In other words, the collision energy is used to vibrationally excite and break the molecule. Higher impact excess energy will lead to a higher kinetic energy, causing the molecule to fragment more easily and to decrease the intensity of molecular ions. We observed that QCEIMS-simulated mass spectra contained fewer fragment ions than their experimental references. For example, the experimental spectrum of 2,4-dimethyl-oxetane (Figure 3a) has 23 product ions, while our QCEIMS simulation produced only four fragment ions plus the molecular ion peak m/z 86. We probed different internal excess energies from 0.2 to 0.8 eV (Figure 4b). With increasing IEE, more fragmentation occurs, increasing the intensity of low mass fragments, but we did not see an increase in the total number of fragments produced. Because the weighted dot-product score gives more weight to the more selective masses found at high m/z ranges, we found that higher IEE values led to decreasing similarity scores. In short, changing ieeatm did not provide a route to improve QCEIMS spectra and we kept the default value of 0.6 eV for subsequent tests.
(3) Initial temperature (tinit)
Last, we investigated the effect of temperature settings ranging from initial temperatures (tinit) of 300 to 1000 K, while keeping all other parameters at default values (Figure 4c). For 2-nonene and adamantane we found that the initial temperatures led to decreasing similarity scores, consistent with the concept that molecules under higher temperature will have more kinetic energy and tend to fragment more easily. For QCEIMS simulations, 2,4-dimethyl-oxetane generated the molecular ion m/z 86 only at low tinit of 300 K, leading to an artificially higher similarity score. As the other two tested molecules also showed their best spectrum similarities at tinit 300K, we chose this parameter value for a final test that utilized a combination of each best setting of ieeatm, ntraj and tinit for each molecule. Interestingly, these simulations did not lead to significant improvements or even to overall decreased similarity scores (see Supporting Information). Therefore, we kept the overall default parameter values for subsequent studies.
Different starting conformers as input for QCEIMS
Local minima on the potential energy surface that are related by rotations around single bonds are called conformational isomers, or conformers. In a mass spectrometer, the conformations of a large cohort of individual chemical molecules are distributed in accord with a Boltzmann distribution at a given ion source temperature. All conformers contribute to the final mass spectrum, to varying degrees related to their relative energies. Ideally, QCEIMS should cover the overall ensemble of conformers. To investigate the impact of the input conformers on the overall QCEIMS results, we selected the highly flexible 2-nonene (PHI=7.51, RBN=5) and the non-flexible adamantane (PHI=1.17, RBN=0) structures. We employed the GMMX software with the Merck Molecular Force Field (MMFF94) to generate starting conformers for individual QCEIMS simulations. For 99 simulations with different starting conformers of 2-nonene, the maximum difference between the lowest-energy and the highest-energy conformer was 2.83 kcal/mol (Figure 5a). For these conformers, dot-product similarity scores ranged from 719 to 824, with a median of 781 and a standard deviation of 24 (Figure 5b). Due to the rigid skeleton and inflexibility of adamantane, GMMX provided only one conformer. Therefore, we used the open source molecular dynamics package CP2K [32] to generate 50 adamantane structures with twisted or stretched bonds that yielded an overall energy range of 5.39 kcal/mol (Figure 5c). Dot-product similarity scores ranged from 849 to 948, with a median similarity of 923 and a standard deviation of 31 (Figure 5d). The examples of these very different molecules showed that QCEIMS similarity scores were independent from input conformer energies (Figure 5a, 5c). Yet, these examples also showed that for both molecules, the QCEIMS fragmentation of specific conformers can lead to quite different dot-product similarities compared to experimental mass spectra, ranging over 100 similarity score units. In addition, we found that dot-product similarities were not normally distributed (Figure 5b, 5d). Our results showed that conformational and other small structural changes may affect QCEIMS simulations. Although adamantane has only a single conformational energy minimum, even slight bond stretches or twists led to quite different mass spectral similarity scores, presumably by biasing molecular dynamics trajectories toward different regions of the potential energy surface. While the QCEIMS software automatically chooses energy-optimized conformers, we propose that a range of different conformers must be calculated to get a good estimate of average mass spectra across the conformational space.
Large scale QCEIMS prediction of small molecule fragmentations
In order to be useful for experimental mass spectrometry, in-silico predictions must not only correctly explain fragmentation and rearrangement reactions for specific molecules, but must also be scalable to generate spectra for hundreds, if not thousands of molecules. Here, we demonstrate the scalability of QCEIMS predictions for small molecules to systematically evaluate parameters and overall accuracies.
The OM2 method only supports carbon, hydrogen, nitrogen, oxygen and fluorine. We therefore chose 451 low molecular weight compounds containing only carbon, hydrogen, nitrogen and oxygen (CHNO). Molecular masses ranged from 26 to 368 Da with an average mass of 129 Da (see Supporting Information). For OM2, computational effort scales as O(N2) ~ O(N3) [33], with N as number of atoms per molecule [33]. The number of single point energy calculations can be estimated to be linearly related to the number of trajectories, and thus linear to the number of atoms. On our computer system with 66 CPU threads, we achieved an average calculation time of 1.55 h per molecule (Figure 6a). Yet, as expected, calculation times exponentially increased with the number of atoms per molecule. For example, with more than 50 atoms, calculation times exceeded 14 hours on the system we had employed (Figure 6a).
Overall, the QCEIMS calculations across all 451 molecules yielded moderately accurate weighted-dot product similarity scores with an average of 608 (Figure 6b) [24]. In GC-MS, similarity scores below 500 are usually not considered for annotation of compounds. While similarity scores above 700 may represent true matches, only scores above 850 are regularly used for direct compound identifications in GC-MS experiments [24]. 47% of all molecules showed good dot-product match factors >700 and 20% of the molecules had excellent scores at >850 similarity. In comparison, lower cosine similarity scores were achieved with an average mass spectral similarity of 557 and a much higher proportion of unacceptably low scoring spectra at similarities <500 (Figure 6b). The regular cosine similarity score does not use weight functions for specific m/z values, unlike the weighted dot product score introduced in 1994 [22] that gives more weight to more specific high m/z product ions in MS fragmentation compared to less specific low m/z fragmentations based on large GC-MS library evaluations. Here, we see a similar trend for QCEIMS spectra.
Molecular descriptors and prediction accuracy
Next we tested the impact of the chemical structures themselves. We used ClassyFire [34] to classify all 451 chemicals into superclasses (supplement file). We found QCEIMS predictions were significantly worse when comparing the organic oxygen superclass of 75 compounds against other superclasses with more than 50 members. Organic oxygen compounds had an average weighted dot-product of 520 whereas the 128 organoheterocyclic compounds achieved significantly better similarities of 648 at p<0.0015 (supplement file). The 100 organic nitrogen compounds yielded an average dot-product similarity of 657 at p<0.001 and the 62 hydrocarbons gave an average of dot-product similarity of 692 at p<0.0001 (supplement file). In conclusion, the QCEIMS method appears to perform worse for oxygen-containing organic compounds than for other major classes. For superclasses with fewer than 50 compounds, statistical tests were deemed to be not robust enough to allow such conclusions.
We also tested if rigid molecules resulted in better prediction accuracy than more flexible ones. Our hypothesis was based on an initial observation that for planar aromatic compounds such as pyridine or aniline, QCEIMS created better quality spectra than for molecules with long chain flexible structures. Our compound data set contained 295 molecules with low flexibility at Kier flexibility index (PHI) < 5 and 161 molecules with high flexibility of PHI > 6. Dot product scores varied significantly across both high-flexibility and low-flexibility molecules (Figure 7a). We found no relationship between flexibility and prediction accuracy. Similarly, we tested rotatable bond number (RBN) as a potential cause for prediction errors (Figure 7b). The median scores for molecules with different RBN values varied between 200-800 and did not depend on increasing RBN. This finding suggests that prediction accuracy is independent of the number of rotatable bonds. In conclusion, we could not find a correlation between flexibility and prediction accuracy at the level of simulation employed.