# Gas-phase reactivity of protonated oxazolone: Chemical dynamics simulations and graph theory-based analysis reveal the importance of ion-molecule complexes

Ariel Perez Mellor, Thomas Bürgi, Riccardo Spezia

► To cite this version:

Ariel Perez Mellor, Thomas Bürgi, Riccardo Spezia. Gas-phase reactivity of protonated oxazolone: Chemical dynamics simulations and graph theory-based analysis reveal the importance of ion-molecule complexes. The Journal of Chemical Physics, 2025, 162 (11), pp.114313. 10.1063/5.0245766. hal-05193419

HAL Id: hal-05193419

https://hal.science/hal-05193419v1

Submitted on 30 Jul 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Gas-phase reactivity of protonated oxazolone: chemical dynamics simulations and graph theory-based analysis reveal the importance of ion-molecule complexes †

Ariel F. Perez Mellor,<fnref n="1" /> Thomas Bürgi,<fnref n="1" /> and Riccardo Spezia<fnref n="2" />

This study delves into the fragmentation mechanisms of the oxazolone form (OXA) of protonated cyclo-di-glycine using chemical dynamics simulations. In this way, we were able to capture the intricate dissociation pathways for different internal energies of the system. By transforming molecular geometries into graph representations, we systematically analyze fragmentation processes and identify key intermediates and ion-molecule complexes (IMCs) that play a crucial role in fragmentation dynamics. The study highlights the distinct isomerization landscapes of OXA, driven by IMC formation, which contrasts with the previously reported behavior of cyclic and linear forms [Perez Mellor et al., J. Chem. Phys., 2021, 155, 124103]. The resulting fragmentation channels are characterized by their distinct energetic thresholds and branching ratios and can provide a molecular explanation of what was observed experimentally. Thanks to an accurate analysis of the trajectories using our graph-theory-based tools, it was possible to point out the particular behavior of OXA fragmentation, which is different from other isomers. In particular, the important role of IMCs is shown, which has an impact in populating different isomeric structures.

## I. INTRODUCTION

Tandem mass spectrometry (MS/MS) is a pivotal analytical technique for examining biological samples due to its high efficiency and standardization.<fnref n="1" />⁻⁴ In proteomics, MS/MS is extensively used, particularly employing collision-induced dissociation (CID) to fragment proteins and peptides into smaller ions.⁵,⁶ The fundamental mechanisms governing peptide fragmentation have been thoroughly investigated through experimental and theoretical studies.³ Resultant peptide fragments are well-characterized, and standardized nomenclatures like the Roepstorff-Fohlman-Biemann system classify these fragments based on cleavage sites and charge localization.⁷,⁸ Among these, singly protonated peptides often generate $b_n^+$ ions, where “b” denotes fragmentation at the peptide bond and n indicates the ion length. These b ions are particularly intriguing because they can cyclize to form macrocyclic structures, potentially leading to internal rearrangements and sequence scrambling, complicating mass spectrometric analyses.⁹⁻¹²

The smallest of the b ions, the $b_2^+$ ion, has attracted significant attention.³,¹³ It has two cyclic forms: the six-membered ring diketopiperazine (DKP) and the five-membered ring oxazolone (OXA). Early experimental and theoretical studies suggested that the oxazolone structure predominates.¹³⁻¹⁶ This was later confirmed using CID coupled with infrared multiple photon dissociation (IRMPD) spectroscopy, which probed the vibrational spectra of $b_2^+$ ions derived from tripeptide fragmentation.¹⁷,¹⁸ However, subsequent studies demonstrated that, for specific peptides, both DKP and OXA structures can be obtained.¹⁹,²⁰

Protonated triglycine, H⁺GGG, is the simplest peptide yielding a $b_2^+$ ion upon fragmentation. Klassen and Kebarle conducted pioneering quantitative studies on the fragmenta-

tion of H⁺GGG.²¹ Siu and colleagues later revisited this system using energy-resolved CID, obtaining threshold energies that aligned well with theoretical predictions.²²,²³ They performed the first quantitative study of $b_2^+$ ion fragmentation.²³ Subsequently, Armentrout and Clark combined threshold CID with *ab initio* calculations of reaction energetics, providing cross sections for the CID of both H⁺DKP and H⁺OXA—the latter produced from *in situ* fragmentation of H⁺GGG. This comprehensive approach revealed intriguing differences between the fragmentation behaviors of H⁺DKP and H⁺OXA ions.²⁴

Chemical dynamics simulations provide a powerful tool for elucidating ion fragmentation mechanisms.²⁵,²⁶ They enable accurate identification of fragmentation products across various species.²⁷,²⁸ This methodology propagates an ensemble of trajectories at energies high enough to observe fragmentation events within feasible computational timescales. Semi-empirical Hamiltonians (e.g., AM1, PM3, PM6-D) or tight-binding density functional theory are typically employed; while they may not offer chemical precision, they provide qualitatively accurate insights.²⁹,³⁰ The analysis of the propagated ensemble based on visual inspection is impractical. Therefore, graph-theory-based methods have been developed and employed. Applications of these methods can be found in protein folding,³¹ enzyme kinetics,³² conformational analysis,³³,³⁴ peptide structure and kinetics studies,³⁵ transition state search,³⁶ and the assignment of vibrational normal modes.³⁷

In this study, we investigate the fragmentation dynamics of protonated oxazolone (H⁺OXA) and compare the results with those of protonated diketopiperazine (H⁺DKP). We apply our recently developed graph theory-based approach to analyze reaction dynamics trajectories.³⁸ This methodology is particularly effective in identifying ion-molecule complexes,

<footnote n="1">¹) Department of Physical Chemistry University of Geneva 30 Quai Ernest-Ansermet, 1211 Geneva 4, Geneva, Switzerland</footnote>

<footnote n="2">²) Laboratoire de Chimie Théorique, Sorbonne Université and CNRS, F-75005 Paris, France</footnote>

<footnote>(*Electronic mail: ariel.perezmellor@unige.ch)</footnote>

<footnote>(Dated: 30 July 2025)</footnote>
elucidating the mechanisms leading to their formation, and understanding how they proceed to form products. Understanding these complexes is crucial for determining roaming mechanisms from a reaction dynamics perspective. Roaming was initially identified by Bowman and co-workers in the photodissociation of formaldehyde³⁹ and has since been observed in numerous reactions.⁴⁰ Recently, roaming was proposed to be responsible for forming certain key products observed in the CID of deprotonated L-cysteine S-sulfate.⁴¹

The article is organized as follows: In Section II, we describe the methodology employed for the chemical dynamics simulations and the techniques developed for result interpretation. Section III presents an overall comparison between H⁺OXA and H⁺DKP, beginning with the sample space of canonical labels (Section III A) and assessing the stability of the graphs based on potential energy distributions extracted from the simulations (Section III B). We compare maximum residence time distributions to probe the formation of ion-molecule complexes (Section III C). An overall comparison using a three-state kinetic model is performed to explain the primary fragmentation of the protonated systems, from which we derive energy thresholds for isomerization and primary fragmentation (Section III D). We then elucidate the isomerization (Section III E) and reaction (Section III F) landscapes of H⁺OXA, highlighting differences compared to H⁺DKP. The article concludes with a summary of our findings in Section IV.

## II. CHEMICAL DYNAMICS SIMULATIONS

In this study, we examine the fragmentation mechanisms of the oxazolone form (OXA⁰⁰) of protonated cyclo-di-glycine (CYC⁰⁰) using chemical dynamics simulations. Figure 1 illustrates the molecular structures. We compare our results with previously reported data on CYC⁰⁰ and its linear intermediate INT⁰⁰³⁸. The superscript "00" designates the activated structure during the chemical dynamics simulation. For simplicity, the protonated moiety (H⁺) has been omitted from the nomenclature in all cases. The simulation protocol mirrors earlier studies, using the RM1-D semi-empirical Hamiltonian⁴² for on-the-fly trajectory propagation. This method, which includes dispersion corrections⁴³, has proven effective in characterizing the dissociation properties of protonated peptides in the gas phase⁴⁴⁻⁴⁶. Fragmentation is induced by introducing excess vibrational energy through microcanonical normal-mode sampling⁴⁷ at various internal energies (127-187 kcal/mol), with rotational energy added via a Boltzmann distribution at 300K. An ensemble of around ten thousand trajectories is propagated for each internal energy. Newton's equations of motion are numerically integrated using the velocity Verlet algorithm and a time step of 0.1 fs, ensuring energy conservation throughout the 40 ps simulations. Table I summarizes the simulation details. The average potential (V) and kinetic (K) energies, along with their standard deviations, were calculated based on the energy decomposition method used in previous studies³⁸. The zero point of the potential energy is set at the electronic energy of the CYC⁰⁰

<figure><img src="image_2.png" /><figcaption>FIG. 1. Schematic representation of the oxazolone form (OXA⁰⁰) of protonated cyclo-di-glycine (CYC⁰⁰). The linear intermediate (INT⁰⁰) formed during the fragmentation of (CYC⁰⁰) is also depicted. The atom numbering used throughout the text is highlighted.</figcaption></figure>

system. All the chemical dynamic simulations are conducted utilizing the software Venus⁴⁸, which seamlessly interfaces with MOPAC 5.022mn⁴⁹ for electronic structure calculations.

<table><caption>TABLE I. Chemical dynamic simulation details. N: Number of trajectories. E<sub>T</sub>: Total energy. E<sub>v</sub>: Activated vibrational energy. ⟨K⟩ and σ⟨K⟩: Mean value and standard deviation of the kinetic energy of the activated structure. ⟨V⟩ and σ⟨V⟩: Mean value and standard deviation of the potential energy of the activated structure.</caption><thead><tr><th rowspan="2">N</th><th rowspan="2">E<sub>T</sub></th><th rowspan="2">E<sub>v</sub></th><th colspan="2">kcal/mol</th></tr><tr><th>⟨K⟩</th><th>σ⟨K⟩</th><th>⟨V⟩</th><th>σ⟨V⟩</th></tr></thead><tbody><tr><td>9999</td><td>138</td><td>127</td><td>63</td><td>10</td><td>75</td><td>10</td></tr><tr><td>9950</td><td>148</td><td>137</td><td>68</td><td>11</td><td>80</td><td>11</td></tr><tr><td>9959</td><td>158</td><td>147</td><td>73</td><td>12</td><td>85</td><td>12</td></tr><tr><td>9945</td><td>168</td><td>157</td><td>78</td><td>13</td><td>90</td><td>13</td></tr><tr><td>9919</td><td>178</td><td>167</td><td>83</td><td>14</td><td>95</td><td>14</td></tr><tr><td>9874</td><td>188</td><td>177</td><td>88</td><td>15</td><td>100</td><td>15</td></tr><tr><td>9785</td><td>198</td><td>187</td><td>93</td><td>15</td><td>105</td><td>15</td></tr></tbody></table>

### A. Trajectory analysis

The ensemble of trajectories was analyzed using a graph theory-based methodology previously developed by some of us.³⁸ Each molecular geometry is represented as an undirected, unweighted, κ-colored graph G = (V, E), where vertices V = {v₁, ..., vₙ} represent atoms, and edges E ⊆ V × V denote potential connections between them. The vertices are color-coded based on atom type, C = {1, 2, ..., κ}. A trajectory is conceptualized as a time-evolving graph, with changes monitored via its adjacency matrix A ∈ R<sup>λ×λ</sup>, where λ represents the number of atoms. This matrix is constructed using threshold bond distances calculated as the sum of atomic radii multiplied by a factor α to account for bond elongation during dynamics. The atomic radii are 0.372Å, 0.781Å, 0.807Å and 0.715Å for hydrogen, carbon, nitrogen, and oxygen, respectively, and α = 1.264.

Threshold bond distances were evaluated by examining the distribution of distances between atom pairs, as shown in Figure S1. This analysis revealed that the first maximum of the distribution falls below the threshold for most atomic pairs, accounting for bond formation and elongation. However, exceptions were noted for specific pairs—hydrogen-hydrogen, nitrogen-nitrogen, nitrogen-oxygen, and oxygen-oxygen—where atom distances exceeded threshold values, indicating the absence of bond formation during the simulation.
tions.

After constructing the adjacency matrix, the canonical label of the graph was computed using the *amtog* and *labelg* tools from the NAUTY package⁵⁰. This process addresses the isomorphism problem, enhancing accuracy in counting and differentiating distinct graphs. Fragments within the graph were identified using the Depth-First Search (DFS) algorithm⁵¹, with additional connectivity and clustering analyses performed based on the threshold distance dclus, which defines clusters as groups of fragments with center-of-mass distances less than 8 Å.

Snapshots were taken every 50 fs, analyzed, and stored in the form of 2D arrays. One such array, CAN[i,j], represents the canonical labels of molecular graphs along trajectory i at discrete time j. The trajectory ensemble analysis focuses on these arrays, and the following properties, explained and discussed in our previous work,³⁸ can be extracted: (1) maximum number of different items, (2) appearance, (3) residence time distribution, (4) maximum residence time distribution, (5) time-dependent probability and standard deviation, (6) mean value or hierarchy, (7) Shannon's entropy, and (8) closure.

## B. New Properties

In addition to the existing properties, two new graph-based properties are introduced here:

9. **Transition Graph:** Consider a directed, weighted, and connected graph $T = (V, E)$, where each vertex in $T$ represents a graph detected during the chemical dynamics simulation. Directed edges between vertices denote transitions between these graphs, and the weights of the edges are given by the transition probabilities $p(i,j)$. The transition probability between two graphs is determined by counting the number of occurrences of the corresponding transition during the simulation normalized by the number of trajectories and snapshots.

It is important to note that constructing $T$ accurately requires capturing simulation snapshots at time intervals shorter than the timescale of processes that induce changes in the adjacency matrix of the system. In our study, proton transfer is the critical process responsible for these changes. Based on this, we selected a quenching interval of 50 fs.

10. **Most Probable Path:** The most probable path $P_{AB}^{mp}$ from vertex A to vertex B is identified by transforming the transition graph $T$ into a graph $T' = (V, E')$, where the edge weights are redefined as $w'(i,j) = [p(i,j) + \epsilon]^{-1}$, with $\epsilon$ being a small positive constant to avoid division by zero. The shortest path between two vertices can be computed using algorithms such as Dijkstra's algorithm⁵².

These two properties are employed to construct and analyze the isomerization and reactivity network. All elements of the sample space $\Omega_{\text{CAN}}$, which correspond to distinct chemical

graphs, form a connected, directed graph. The transition probabilities between these graphs are derived from the chemical dynamics simulation, providing insight into the isomerization pathways and chemical reactivity.

## C. Statistical analysis and kinetic models

The primary fragmentation of protonated cyclodipeptides is modeled using a three-state framework, where each snapshot from the chemical dynamics simulation is categorized into one of three states: the starting point (SP), intermediate (INT), and primary fragmentation (PF) states. The SP state consists of graphs resembling the activated structure. In contrast, the PF state includes disconnected graphs that remain so for a specific duration $\Delta$, ensuring they do not re-enter the connected region. The INT state encompasses all other snapshots. The population of these states over time can be analyzed using either full simultaneous fitting (called here and hereafter Method A) or constrained simultaneous fitting (Method B), with the latter focused on determining the most relevant reactive rate constant via direct flux measurements.⁵³ Further details on these methodologies are available in our previous work.³⁸ From this analysis, we obtained microcanonical rate constants, and then threshold energies are calculated from Rice-Ramsperger-Kassel-Marcus (RRKM) theory⁵⁴. This approach facilitates comparisons of dynamic behaviors across different systems.

## III. RESULTS AND DISCUSSION

The results for the OXA⁰⁰ system are presented and discussed, followed by a comparative analysis with the CYC⁰⁰ and LIN⁰⁰ systems. For CYC⁰⁰ and LIN⁰⁰, the data were obtained from a previous study³⁸ where the same methodology for chemical dynamics simulation was applied.

### A. The sample space of CAN array

The analysis begins by examining the sample space of the CAN array ($\Omega_{\text{CAN}}$), which consists of all graphs visited at a given total energy and for a specified number of trajectories. We have chosen the same total energy of 178 kcal/mol to ensure a meaningful comparison between the systems. The sample space $\Omega_{\text{CAN}}$ depends on the total energy (E), the number of trajectories (N), and the trajectory length ($\Lambda$)), and can be divided into connected and disconnected regions. Figure 2 illustrates a schematic representation of the sample space for the oxazolone system. The connected region corresponds to where isomerization processes occur, while the disconnected region is associated with fragmentation and ion-molecule complex formation. The graphs are labeled using the permutation-invariant canonical label in Graph6 format. Additionally, the characters S, I, and P were prefixed to the labels for further discrimination, as will be detailed later. It should be noted that some graphs are not chemically possible
<figure><img src="image_3.png" /><figcaption>FIG. 2. Schematic representation of the sample space ($\Omega_{CAN}$) of OXA$^{00}$, which is divided into two regions: connected and disconnected. The graphs are labeled using the canonical label in the Graph6 format. Hydrogen atoms are highlighted in green, carbon atoms in gray, nitrogen in blue, and oxygen in red. The sample space has been partitioned using a three-state model (see text for further details), with a color code denoting the states: light blue represents the starting point, light gray corresponds to the intermediate state, and light red indicates the primary fragmentation state.</figcaption></figure>

and have been intentionally highlighted with a question mark symbol. They also belong to the sample space.

Figure 3, panels a-c, shows the average time-dependent probability of each canonical label (plotted as index) within the sample space, distinguishing between connected (red dots) and disconnected (blue dots) graphs. For the OXA$^{00}$ system, the sample space includes 602 graphs, with disconnected graphs being more populated. The most populated connected graph is highlighted with a blue arrow in the inset of panel a. The CYC$^{00}$ and LIN$^{00}$ systems display significantly more graphs (1498 and 1372, respectively), with connected graphs being the most frequent in these systems. Panels d-f of Figure 3 depict Shannon entropy values for CAN arrays as a function of the number of trajectories, showing stabilization after about 2000 trajectories. Shannon entropy quantifies the uncertainty or diversity in the distribution of canonical labels (graphs) within the sample space and is defined as $S_{CAN} = -\sum_i p_i \ln(p_i)$, where $p_i$ is the mean value of the time-dependent probability of occurrence of the i-th canonical label. A lower Shannon entropy value indicates that the system's behavior is dominated by a few specific graphs or pathways, reflecting a more ordered or predictable fragmentation process. In contrast, a higher Shannon entropy value signifies a more uniform distribution among many different graphs, indicating a greater diversity of isomerization/fragmentation pathways and a higher degree of disorder in the system.

In this context, the lower entropy of the OXA$^{00}$ system (1.77) suggests that it predominantly follows limited isomerization/fragmentation pathways. The higher entropy values for the CYC$^{00}$ and LIN$^{00}$ systems (2.73 and 2.95, respectively) reflect a broader range of accessible pathways and greater structural diversity in their isomerization/fragmentation behavior. The similarities in CYC$^{00}$ and LIN$^{00}$ systems are attributed to a common mechanistic path-

<figure><img src="image_4.png" /><figcaption>FIG. 3. Left Panel: Analysis of the average time-dependent probability for each canonical label plotted as index in (a) OXA$^{00}$, (b) CYC$^{00}$, and (c) LIN$^{00}$ at a total energy of 178 kcal/mol. The connected graphs are in blue, whereas the disconnected ones are in red. The insets provide a zoomed view of the first 40 canonical labels. Right Panel: Evaluation of the Shannon entropy of the canonical label array as the function of the number of trajectories in (d) OXA$^{00}$, (e) CYC$^{00}$, and (f) LIN$^{00}$.</figcaption></figure>

way involving initial ring-opening and proton transfer,$^{38}$ leading to similar behavior between the two systems.

## B. Potential energy distribution and rotational constants of a given graph

Each canonical label is characterized by various attributes, including kinetic and potential energy distributions, rotational constants, and shape descriptors, among others. These attributes are essential for the identification of conformers associated with a specific molecular graph and for assessing its stability. Figure 4, panel a-c, presents the potential energy distribution for the three systems under the same total energy. The CYC$^{00}$ graph, the most stable configuration, has a distribution centered at 90.12 kcal/mol, with the potential energy zero set at the electronic energy minimum of CYC$^{00}$. The OXA$^{00}$ graph, less stable by 4.6 kcal/mol, exhibits a higher
<figure><img src="image_5.png" /><figcaption>FIG. 4. Left Panel: Potential energy distribution of (a) OXA<sup>00</sup>, (b) CYC<sup>00</sup>, and (c) LIN<sup>00</sup> at the same total energy of 178 kcal/mol. Right Panel: Rotational constants distribution for (d) OXA<sup>00</sup>, (e) CYC<sup>00</sup>, and (f) LIN<sup>00</sup> at a total energy of 178 kcal/mol, sorted in ascending order as A > B > C. At t = 0, only the distribution of the largest rotational constant is shown.</figcaption></figure>

energy distribution. The least stable system, LIN<sup>00</sup>, is 11.5 kcal/mol higher than CYC<sup>00</sup>. The data presents an intriguing, non-linear relationship between the Shannon entropy of the CAN array and their stability. The LIN<sup>00</sup> exhibits the highest Shannon entropy among the studied graphs. Furthermore, the CYC<sup>00</sup>, characterized by its substantial Shannon entropy, ranks as the most stable graph. This stability surpasses that of OXA<sup>00</sup>, which, despite being less energetically stable, has a lower Shannon entropy than OXA<sup>00</sup>.

The number of conformers for each graph is inferred from the distribution of rotational constants. For each canonical label, the rotational constants are calculated from the eigenvalues of the moments of inertia.⁵⁵ Figure 4, right panel, displays the rotation constant distribution for t = 0 and t > 0, sorted in ascending order (A>B>C). Only the distribution of A is shown for t = 0. In the OXA<sup>00</sup> system, a single conformer is selectively activated at t = 0, with a distribution centered at 0.18 cm<sup>-1</sup>. During the chemical dynamic simulations, another conformer activates, distinguished by a distribution centered at 0.14 cm<sup>-1</sup>, corresponding to free rota-

<figure><img src="image_6.png" /><figcaption>FIG. 5. Projection on a plane of atoms' position for OXA<sup>00</sup> (top), CYC<sup>00</sup> (middle), and LIN<sup>00</sup> (bottom) at E<sub>T</sub> = 178kcal/mol. The left panel represents t = 0, while the right encompasses t > 0. Hydrogen atoms are highlighted in green, carbon atoms in gray, nitrogen in blue, and oxygen in red. The colorbar representing occurrences is qualitatively set to indicate low (L) in blue, medium (M) in orange, and high (H) in pink.</figcaption></figure>

tion around the C<sub>3</sub>–C<sub>4</sub> bond. The geometrical oscillations between these conformers correspond to two positions of the dihedral angle N<sub>2</sub>C<sub>3</sub>C<sub>4</sub>N<sub>5</sub>, specifically at 0° (cis) and 180° (trans). Figure 5 (top) shows the projection of the position of the atoms on a plane formed by atoms C<sub>1</sub>, N<sub>2</sub>, and O<sub>6</sub> using TRAVIS software.⁵⁶ The CYC<sup>00</sup> system shows a similar distribution between t = 0 and t > 0, centered near 0.15cm<sup>-1</sup>, indicating the activation of either a single conformer or two highly similar conformers. The latter scenario is indeed the case. These conformers are differentiated by the orientation of the hydroxy group, which alternates between cis and trans configurations. Figure 5 (middle) highlights this observation. This motion induces minimal changes in the moments of inertia of the molecule and, consequently, has a negligible effect on the rotational constant distribution. The LIN<sup>00</sup> system has a single conformer at t = 0, characterized by a rotational constant distribution centered at 0.22 cm<sup>-1</sup>. As chemical dynamic simulations progress, long-amplitude movements activate a new conformer, with a frequency distribution centered at 0.11cm<sup>-1</sup>. Figure 5 (bottom) summarizes this observation.
### C. Maximum residence time and clustering

Figure 6, panels a-c, shows the maximum residence time (MRT) distribution for the primary fragmentation guess (PFG) matrix, MRT<sub>PFG</sub>. The MRT denotes the duration for which the graphs persist continuously within the disconnected region without returning to the connected domain. This metric is a valuable estimate for the Δ parameter. The latter corresponds to the time we have to look ahead to classify a given snapshot as a primary fragmentation. The probability density function (PDF) characterizing the MRT exhibits similar trends for the CYC<sup>00</sup> and LIN<sup>00</sup> systems, where a monotonically fast decay is observed. On the other hand, the OXA<sup>00</sup> system shows a more complex behavior. There is a rapid drop up to approximately 1 ps. This phase is followed by a notable increase, reaching a peak at 1.55 ps (red dashed line). Subsequently, the PDF exhibits an exponential decay. Such behavior suggests the formation of long-lived ion-molecule complexes (IMCs) with the mass of the precursor ion, as we will see later. In this case, we have selected Δ = 15 ps (black solid line) to ensure that only 0.01% of the trajectories are incorrectly assigned to the PF state. This Δ value is three times larger than that employed for CYC<sup>00</sup> and LIN<sup>00</sup> systems (black dashed lines) where the formation of IMCs (with the mass of the precursor ion) is negligible.<sup>38</sup>

In the present case, however, the formation of IMC is obtained, as visible from pair distribution distances between fragment centers of mass. Notably, Figure 6 (panels d-f), shows center-of-mass distances among fragments in OXA<sup>00</sup> (d), CYC<sup>00</sup> (e), and LIN<sup>00</sup> (f). The fragments were detected by applying the DFS algorithm to the adjacency matrix, and their mass was calculated by summing the atomic masses of all atoms belonging to each fragment. The figure employs a dual-layered graphical representation. The underlying layer, shaded in gray, illustrates the distribution encompassing all mass fragments without any specific categorization. Overlaid on this is the breakdown into principal components in a palette of green, cyan, blue, and red. In the study of the oxazolone system, two main contributions are identified, represented by a set of masses of {30.03, 29.03, 27.99, 27.99} and {85.02, 30.03}. The first contribution exhibits a peak at 3.58Å with a full width at half maximum (FWHM) of 1.33Å. This peak is indicative of the formation of an IMC with m/z of 59.06 involving the methaniminium cation CH<sub>2</sub>NH<sub>2</sub><sup>+</sup> (m/z 30.03) and the neutral methylene imine CH<sub>2</sub>NH (m/z 29.03), which occurs after the release of two carbon monoxide molecules from the parent OXA<sup>00</sup>. This IMC formation was experimentally proposed and validated<sup>23,24</sup>. Conversely, the second contribution is delineated by the emergence of a different IMC involving the same methaniminium cation and the residual five-member ring C<sub>3</sub>H<sub>3</sub>NO<sub>2</sub> (m/z 85.02). The m/z of the complex is identical to that of OXA<sup>00</sup>, 115.05. This IMC formation is marked by a maximum at 4.34 Å and a FWHM of 1.48 Å. In this scenario, the IMC materializes directly due to the cleavage of the C<sub>3</sub>–C<sub>4</sub> bond of the parent molecule, leading to a five-member ring (oxazolone structure) incorporating a carbene entity. The formation of this IMC is discussed in more detail later in this article.

<figure><img src="image_7.png" /><figcaption>FIG. 6. Left Panel: Distribution of the maximum residence time in the primary fragmentation guess matrices of (a) OXA<sup>00</sup>, (b) CYC<sup>00</sup>, and (c) LIN<sup>00</sup> at the same total energy of 178 kcal/mol. The line representing the chosen delta value is distinctly emphasized for each system. Right Panel: Probability density function of the distance between disconnected graphs in (d) OXA<sup>00</sup>, (e) CYC<sup>00</sup>, and (f) LIN<sup>00</sup> at a total energy of 178 kcal/mol. The graph features a black line indicating the designated dclus at 8 Å. The distribution is normalized to its maximum value for ease of comparison. Furthermore, the breakdown into principal components is delineated for each system.</figcaption></figure>

The CYC<sup>00</sup> and LIN<sup>00</sup> systems share the same main contribution to IMC formation. They are {58.03, 29.03, 27.99}, {57.02, 30.03, 27.99} and {30.03, 29.03, 27.99, 27.99}. The latest has the same nature as in OXA<sup>00</sup>, with a maximum at 3.66 Å. The first two correspond to an IMC formed after the release of a CO molecule (m/z 27.99). In this case, the methaniminium cation and the methylene imine are bound to the linear C<sub>2</sub>H<sub>3</sub>NO and C<sub>2</sub>H<sub>4</sub>NO<sup>+</sup>. In each scenario, the peak of the distribution curves consistently falls below the predetermined cut-off cluster distances (dclus = 8 Å). This alignment is essential as it confirms the accuracy of cluster formation estimates, which are critical for determining the reaction mechanisms.
## D. Three-State Kinetic Model in OXA<sup>00</sup>

In this section, we analyze the fragmentation kinetics of OXA<sup>00</sup> using the previously defined Δ parameter. The model considers three distinct states: starting point (SP), intermediate (INT), and primary fragmentation (PF). The SP state includes graphs that resemble the activated structure, whereas the PF state consists of disconnected graphs that remain disconnected for a duration Δ, indicating they do not re-enter the connected region (i.e., primary fragmentation products). The INT state represents all remaining graphs. Figure 2 summarizes the partitioning of the sample space of the canonical labels array according to the three-state model. The same canonical label may be assigned to two different cases, as is the case with certain IMC graphs, depending on whether they are bound or unbound. To address this issue and facilitate tracking, we prefix the canonical labels with the characters S (starting point), I (intermediate), and P (primary fragmentation). This approach allows us to distinguish between graphs that form IMCs in the mass of the parent channel and those that undergo primary fragmentation.

Figure S2 shows graphs that encompass the SP state of the OXA system, with the starting activated graph (SAG) being the primary component CAN<sup>1</sup>. The probabilities of connected and disconnected corrections are significantly lower than those of SAG, indicating a minimal contribution to the SP probability. The time evolution of the probability of SP, INT, and PF states obtained from the chemical dynamic simulation at different vibrational energies is depicted in Figure 7 in gray-colored lines. The SP state decreases rapidly, the INT state increases to a maximum before decreasing, and the PF state steadily increases.

To quantify the evolution of the states and extract effective rate constants, we employ a three-state kinetic model, see the bottom left panel in Figure 7, governed by the following differential equations:

$$ \frac{dN_{SP}}{dt} = k_{SI}N_{INT} - (k_{IS} + k_{PS})N_{SP} $$

$$ \frac{dN_{INT}}{dt} = k_{IS}N_{SP} - (k_{SI} + k_{PI})N_{INT} $$

$$ \frac{dN_{PF}}{dt} = k_{PS}N_{SP} + k_{PI}N_{INT} $$

Here, N<sub>SP</sub>, N<sub>INT</sub>, and N<sub>PF</sub> represent the populations of the SP, INT, and PF states, respectively. The rate constants k<sub>IS</sub>, k<sub>SI</sub>, k<sub>PS</sub>, and k<sub>PI</sub> govern the transitions between these states:

* k<sub>IS</sub>: Rate constant for the transition from SP to INT, representing the effective isomerization step.

* k<sub>SI</sub>: Rate constant for the reverse transition from INT to SP, accounting for reversion to the initial state.

* k<sub>PS</sub>: Rate constant for direct fragmentation from SP to PF, bypassing the intermediate state.

* k<sub>PI</sub>: Rate constant for fragmentation from INT to PF.

By solving these coupled differential equations, the time-dependent populations of the three states can be determined.

<figure><img src="image_8.png" /><figcaption>FIG. 7. Populations of the three states (represented by gray lines) at a given vibrational energies (measured in kcal/mol) for the OXA<sup>00</sup> system. The colored lines correspond to the fitted functions (Method A), with blue representing the starting-point state, black representing the intermediate state, and red representing the primary fragmentation state. The bottom left panel highlights the kinetic scheme followed for the simultaneous fitting.</figcaption></figure>

The solution is obtained using the Laplace-Carson transform method.<sup>57</sup> For extracting the effective rate constants, we employed two approaches:<sup>38</sup>

1. **Simultaneous Fit (Method A):** All rate constants are obtained from the simultaneous fitting using the solution of the three-state model. There are three populations to be fitted and four parameters to be extracted.

2. **Flux-Based Fit (Method B):** The primary fragmentation rate constant (k<sub>PI</sub>) is obtained from the reaction flux<sup>53</sup> and its value is kept fixed during the simultaneous fitting. Figure S3 summarizes the extraction of k<sub>PI</sub>, defined as the value at which the maximum in the dis-
<figure><img src="image_9.png" /><figcaption>FIG. 8. Microcanonical rate constants (symbols) for various processes plotted against their corresponding internal energies for OXA<sup>00</sup>. The colored line illustrates the fit derived from the RRKM model, where parameters such as energy threshold and effective frequency factor are extracted. Panel (a) shows results obtained from the simultaneous fitting method, while panel (b) shows those from the flux method.</figcaption></figure>

tribution is observed (red line). In this case, we have three populations to be fitted and three parameters to be extracted.

Both methods show a strong correlation with the observed evolution of the three states, accurately capturing their dynamics over time. Refer to the colored lines in Figure 7 for Method A and in Figure S4 for Method B. The rate constants obtained from both methods are summarized in Table S1.

Figure 8 shows the evolution of the rate constants as a function of vibrational energy. Both methods yield similar trends for all rate constants, except for $k_{SI}$. Method A results in an increase in $k_{SI}$ as a function of energy, while Method B shows minimal variation. The difference in $k_{SI}$ can be attributed to the interdependence between $k_{SI}$ and $k_{IS}$. Interdependence refers to the relationship between these two rate constants, where changes in one can influence the other. For Method A, there may be a small leakage effect, meaning that some of the influence from $k_{IS}$ might inadvertently affect the estimation of $k_{SI}$, leading to less accurate results. In contrast, in Method B, by fixing one of the rate constants, the fitting is constrained to three parameters, thereby reducing the potential for leakage and improving the reliability of the fitted values.

The energy dependence of the rates follows the RRK theory<sup>54</sup>, which allows the determination of activation energy by fitting the following equation:

$$k(E_v) = v \left(1 - \frac{E_a}{E_v}\right)^{s-1} \quad (1)$$

where $E_a$ is the activation energy, $E_v$ is the vibrational energy, $s$ is the number of vibrational degrees of freedom, and $v$ is a frequency factor. As shown in Figure 8, this expression fits the data well, consistent with previous studies<sup>30,38,58</sup>. The threshold energies obtained are summarized in Table S2. The energy threshold for isomerization is approximately 11 kcal/mol, while the threshold for direct fragmentation is about 21 kcal/mol. The energy required for the reaction from an effective intermediate is 15 kcal/mol. The relatively low energy threshold for isomerization is attributed to the formation of the

<figure><img src="image_10.png" /><figcaption>FIG. 9. Comparison of effective energy thresholds between systems: CYC<sup>00</sup>, LIN<sup>00</sup>, and OXA<sup>00</sup>.</figcaption></figure>

IMC, primarily involving the cleavage of the $C_3-C_4$ bond. This point will be discussed further in a later section.

### Comparison with the CYC<sup>00</sup> and LIN<sup>00</sup> systems

The use of the three-state model to quantify the fragmentation of CYC<sup>00</sup> and LIN<sup>00</sup> has been previously investigated<sup>38</sup> under the same theoretical approximations as those applied to the OXA<sup>00</sup> system. In the CYC<sup>00</sup> system, LIN<sup>00</sup> emerges as the most populated isomer, formed through ring opening via cleavage of the $C_1-C_2$ bond, followed by proton transfer from $N_6H$ to $C_2$ (see Figure 1 for the atom numbering). Consequently, these two systems share a largely similar fragmentation basin. At a total energy of 178 kcal/mol, the PF state exhibits a potential energy distribution centered at 118.02 kcal/mol for CYC<sup>00</sup> and 119.24 kcal/mol for LIN<sup>00</sup>, see Figure S5. The OXA<sup>00</sup> system, however, has a potential energy distribution centered at 115 kcal/mol, which is slightly different. For the SP state, CYC<sup>00</sup> has the lowest energy, centered at 90.21 kcal/mol, while LIN<sup>00</sup> is centered at 101.64 kcal/mol. The OXA<sup>00</sup> system's SP state is 6 kcal/mol lower than LIN<sup>00</sup> and 5 kcal/mol higher than CYC<sup>00</sup>. According to the three-state model, the OXA<sup>00</sup> system exhibits the lowest barrier for isomerization at 10.6 kcal/mol, compared to 30.5 kcal/mol for CYC<sup>00</sup> and 13.8 kcal/mol for LIN<sup>00</sup>, as seen in Figure 9. This explains why the OXA<sup>00</sup> system displays the highest isomerization rate constant, $182.4 \cdot 10^{-5} \text{fs}^{-1}$ compared to $8.6 \cdot 10^{-5} \text{fs}^{-1}$ for CYC<sup>00</sup> and $13.4 \cdot 10^{-5} \text{fs}^{-1}$ for LIN<sup>00</sup>. Regarding the energy required to react from an effective intermediate state, the OXA<sup>00</sup> system requires only 15.5 kcal/mol, whereas CYC<sup>00</sup> and LIN<sup>00</sup> require 37.9 kcal/mol and 26.4 kcal/mol, respectively. The difference in energy between OXA<sup>00</sup> and CYC<sup>00</sup> is approximately 0.9 eV, which is comparable to the energy threshold reported experimentally, 1.0 eV.<sup>24</sup> Another difference is that the $k_{PS}$ rate constant is significantly higher in the OXA<sup>00</sup> system than in the CYC<sup>00</sup> system. In the former, the reaction from the starting point proceeds via the cleavage of the $C_3-C_4$ bond, leading to the direct formation of the methaniminium cation.

Although the energy threshold obtained here cannot be directly compared to experimental values due to the use of a large excess of vibrational energy to induce fragmentation events on the picosecond timescale, it provides a qualitative estimation when comparing different systems, as done in this study.
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 10. a. Isomerization network of OXA<sup>00</sup>, comprising the eight most populated connected structures at a total energy of 178 kcal/mol. The radii of the nodes in this network are proportionally scaled to reflect the self-transition probability of each structure, and the width of the arrows is proportional to the transition probability. The color code depicts the initial state in blue and the intermediate states in gray. The values displayed represent the average potential energy measured in kcal/mol. The four tiny dots denote the formation of IMC. b. Most probable path that connect OXA<sup>00</sup> to LIN<sup>02</sup>. The rectangular red box highlights the formation of IMCs.</figcaption></figure>

E. Isomerization landscape in OXA<sup>00</sup>

The isomerization landscape of the OXA<sup>00</sup> system was investigated using the graph theory methodology. The landscape is divided into two regions: the connected region, where the isomers are found, and the disconnected region, where IMC formation occurs. The eight most populated connected structures are depicted in Figure 10a. We employed a representation that captures two key aspects: first, the isomer abundance, which is computed through the self-transition probability of each graph, and second, the transition probability between them. They were determined using the transition graph property, which is one of the new properties presented in this article. The radii of the nodes are proportional to the isomer abundance, while the width of the arrows is proportional to the transition probability. The color code depicts the initial state in blue and the intermediate states in gray. When a

<figure><img src="image_13.png" /><figcaption>FIG. 11. Isomerization network emphasizing the formation and inter-connectivity of the most significant ion-molecule complexes at a total energy of 178 kcal/mol for OXA<sup>00</sup>. The values shown on the network correspond to the mean potential energy, quantified in kcal/mol. The four tiny dots indicate the formation of IMC.</figcaption></figure>

chemically feasible structure exists, it is represented by a 2D chemical formula; otherwise, a graph representation is used. The connected isomers are further classified into three categories: cyclic isomers, such as CYC<sup>01</sup>, CYC<sup>02</sup>, and CYC<sup>03</sup>; linear structures, such as LIN<sup>01</sup> and LIN<sup>02</sup>; and turning points, TUR<sup>01</sup> and TUR<sup>02</sup>. The latter will be discussed in more detail subsequently.

The most abundant connected isomer is LIN<sup>01</sup>, which corresponds to ring opening through the cleavage of the C<sub>7</sub>–C<sub>6</sub> bond (see Figure 1 for atom numbering). Surprisingly, the second most populated connected graph is TUR<sup>01</sup>, which is not chemically feasible due to a hydrogen atom being bonded to both carbon and nitrogen atoms (this is underlined by a red question mark in the figure). As we will emphasize later, this type of graph is associated with the formation of IMC. Energetically, it corresponds to a situation in which the total energy is equal to the potential energy, and thus we called these structures, here and hereafter, "turning points". Finally, the most abundant cyclic structure, apart from the activated OXA<sup>00</sup>, is CYC<sup>01</sup>, which involves a proton transfer from one of the protons at C<sub>4</sub> to C<sub>3</sub>, making the C<sub>3</sub> carbon atom asymmetric. Intriguingly, no pathway within the connected region links the starting graph with the second most populated linear structure, LIN<sup>02</sup>. Based on chemical intuition, determining the possible mechanism is exceedingly challenging. Consequently, we used the "most probable path" property, as illustrated in Figure 10b. The mechanism involves the formation of IMCs highlighted with the rectangular red box. Initially, there is cleavage of the C<sub>3</sub>–C<sub>4</sub> bond forming IMC, followed by a proton transfer that leads to the formation of a carbocation. Subsequently, ring opening occurs, and the neutral methylene imine reattaches, forming the LIN<sup>02</sup> structure.

Our simulations reveal that IMCs play a crucial role in the isomerization mechanism during the fragmentation of OXA<sup>00</sup>.
<figure><img src="image_14.png" /><figcaption>FIG. 12. Plane-projection probabilities (left panel) and the distribution of distances between the center of mass (right panel) in the ion-molecule complex IMC<sup>01</sup> for OXA<sup>00</sup> at a total energy of 178 kcal/mol. The colored regions depict the probability of locating the projection of the residues, NH<sub>2</sub>–CH<sub>2</sub><sup>−</sup> or NH–CH<sub>2</sub>, at specific points on a plane defined by the C<sub>1</sub>N<sub>2</sub>O<sub>6</sub> atoms in the structure. Hydrogen atoms are highlighted in green, carbon atoms in gray, nitrogen in blue, and oxygen in red. The colorbar representing occurrences is qualitatively set to indicate low (L) in blue, medium (M) in orange, and high (H) in pink.</figcaption></figure>

Figure 11 illustrates the isomerization network, which is composed exclusively of IMCs. The most prevalent isomer within this network is IMC<sup>01</sup>, which exhibits a mean potential energy of 102.25 kcal/mol. In this configuration, the methaniminium cation exhibits a dynamic "roaming" behavior around the five-membered ring fragment, eventually initiating a proton transfer event. This transfer enables the formation of additional IMC isomers, such as IMC<sup>01</sup>. The turning points in the isomerization network, represented as TUR<sup>01</sup>, TUR<sup>02</sup>, and TUR<sup>03</sup>, are identified just before the proton transfer event. The proton is shared between the two molecular moieties. These turning points mark the onset of proton migration.

In this study, IMCs are formed through charge-dipole interactions. To illustrate this phenomenon, we have plotted the positions of either the methaniminium cation or the neutral methylene imine across all trajectories, projected onto the plane defined by the cyclic structure, as shown in Figure 12a-

<figure><img src="image_15.png" /><figcaption>FIG. 13. Branching ratios at a total energy of 178 kcal/mol for OXA<sup>00</sup>, with the left panel depicting the breakdown of the intermediate into various classes: linear, cyclic, and ion-molecule complex structures. The right panel focuses on the major contributions to primary classes, specifically the ion-molecule complex and linear structures.</figcaption></figure>

c. Snapshots corresponding to the same canonical label were selected, ensuring that none of the structures met the criteria for chemical reactivity. This criterion was assessed by focusing on canonical labels indexed with an "I" for intermediate states (see Figure 2). For the IMC<sup>01</sup> isomer, shown in panel a, the methaniminium cation exhibits a large-amplitude motion around the remaining five-membered ring structure. The amplitude of this motion can exceed 10 Å, with the maximum displacement localized at 4.38 Å (see panel d of the same figure). This roaming motion is highly asymmetric and tends to localize around specific atoms within the cyclic structure, particularly C<sub>3</sub>, O<sub>6</sub>, C<sub>7</sub>, and O<sub>8</sub>. These asymmetrical interactions suggest that the charge-dipole interaction driving the motion is not uniformly distributed across the structure, with the identified atoms acting as proton acceptors, stabilizing the roaming cation. In the case of the IMC<sup>01</sup> complex, see panel b, the neutral methylene imine moves around the five-membered ring, which remains protonated. The motion is more localized around the protonated site, particularly at C<sub>3</sub>. The center-of-mass distance distribution shows a peak at a shorter distance of 4.49 Å, panel e, with a significantly narrower width compared to IMC<sup>01</sup>, indicating more confined movement in this complex. For the IMC<sup>01</sup> complex, panel c, the methaniminium cation roams around an isomer of the five-membered ring, where the charge is stabilized by proton transfer from N<sub>2</sub>. The amplitude of the motion is centered at 4.61 Å, panel f, with a distribution width of 1.10 Å, falling between the values observed for IMC<sup>01</sup> and IMC<sup>02</sup>. In this case, the movement is primarily localized around N<sub>2</sub>, which acts as the proton acceptor. The distribution width of the center-of-mass distance provides insight into the degree of localization of the motion during IMC formation, with narrower distributions indicating more constrained motion within the complex.

Finally, the role of IMCs in the chemical dynamics simulation is further explored by decomposing the population of the intermediate state into linear, cyclic, and IMC structures and calculating their relative populations over time with respect to the total population of the intermediate state. Figure 13a provides a summary of these results. During the first 2 ps, the
<figure><img src="image_16.png" /><figcaption>FIG. 14. Principal isomerization network of CYC$^{00}$, comprising the eight most populated connected structures at a total energy of 178 kcal/mol. The radii of the nodes in this network are proportionally scaled to reflect the self-transition probability of each structure. The values displayed represent the average potential energy measured in kcal/mol.</figcaption></figure>

intermediate state population is primarily composed by IMCs and linear structures. As the simulation progresses, the IMC structures begin to dominate, eventually exceeding 90% of the total intermediate population. In contrast, cyclic isomers, such as $CYC_{01}^{01}$, remain comparatively less prevalent, reaching a maximum of approximately 20% at the end of the dynamics. Panel b of Figure 13 illustrates that the $IMC_{01}^{01}$ structure is the most frequently occurring ion-molecule complex, while $LIN^{01}$ is the most populated linear structure. These two structures dominate the isomerization landscape throughout the chemical dynamics simulation in $OXA^{00}$ system.

### Comparison with the CYC$^{00}$ and LIN$^{00}$ systems

The isomerization landscapes of the CYC$^{00}$ and LIN$^{00}$ systems are largely similar. Both systems share the same eight most populated graphs, differing only in their population (see Figures 14 and 15). The mean values of the potential energy distributions are statistically consistent between the datasets, with identical means within the margin of error. In both cases, linear structures dominate the networks, driven by the proton transfer dynamics after ring opening. As previously discussed,$^{38}$ the primary isomerization mechanism involves ring opening due to cleavage of the $C_1-C_2$ bond, followed by proton transfer. It is important to note that all linear structures in the networks represent tautomeric forms.

A significant observation is the absence of turning points in both CYC$^{00}$ and LIN$^{00}$, which suggests that the formation of IMCs in the parent ion mass is highly improbable in these

<figure><img src="image_17.png" /><figcaption>FIG. 15. Principal isomerization network of LIN$^{00}$, comprising the eight most populated connected structures at a total energy of 178 kcal/mol. The radii of the nodes in this network are proportionally scaled to reflect the self-transition probability of each structure. The values displayed represent the average potential energy measured in kcal/mol.</figcaption></figure>

systems, as discussed earlier.$^{38}$ This contrasts with the $OXA^{00}$ system, where IMC formation is a key feature of the isomerization landscape. The major distinction between the isomerization networks of $OXA^{00}$ and $CYC^{00}$ is the predominance of IMC formation in $OXA^{00}$, while in the $CYC^{00}$ and $LIN^{00}$ systems, the occurrence of IMCs is minimal.

One structure, however, is common between the isomerization networks of $OXA^{00}$ and $CYC^{00}$: the $LIN^{02}$ isomer in $OXA^{00}$ corresponds to $LIN^{04}$ in the $CYC^{00}$ isomerization network. In the $OXA^{00}$ system, we have demonstrated that $LIN^{02}$ is populated via IMC formation, as shown in Figure 11b. This observation reveals the isomerization pathway between oxazolone and diketopiperazine structures, and vice versa, involving the formation of IMC intermediates.

## F. Fragmentation Mechanisms

Fragmentation mechanisms are analyzed by assessing branching ratios among mass distributions that meet fragmentation criteria and the probability of the primary fragmentation state. This approach allows the categorization of the fragmentation process into four distinct types: primary fragmentation, secondary fragmentation, tertiary fragmentation, and non-assigned events. Primary, secondary, and tertiary fragmentation refers to the dissociation of the molecule into two, three, and four clusters, respectively. The non-assigned events encompass cases where more than two fragments are formed, but the inter-fragment distances are lower than 8 Å, rendering them challenging to categorize.
<figure><img src="image_18.png" /><figcaption>FIG. 16. Branching ratios for OXA⁰⁰ at a total energy of 178 kcal/mol. (a) Decomposition of the primary fragmentation state into distinct reaction channels. (b) Relative population of the *m/z* channels. (c) Analysis of the CO loss channel, delineating direct and ion-molecule formation pathways. (d) Analysis of the neutral methylene imine loss channel, distinguishing between direct and ion-molecule formation pathways.</figcaption></figure>

For OXA⁰⁰ system, two main primary fragmentation pathways are observed. The first pathway leads to the formation of the *m/z* 87.06 (C₃H₇N₂O⁺) fragment, resulting from the loss of a CO molecule. By analyzing its formation as a function of time (see Figure 16a), we notice that this fragment is obtained in a very short time-scale (a maximum relative abundance of about 60% is reached within the first picosecond), and then it decreases gradually down to about 1%. The second primary fragment, the methaniminium cation (m/z 30.03), is obtained via the loss of a neutral C₃H₃NO₂. Initially, its abundance decreases to around 25% before stabilizing at approximately 30%. A minor contribution, about 4%, to the primary fragmentation process is also observed, corresponding to the formation of *m/z* 88.04 (C₃H₄NO₂⁺) after the release of neutral methylene imine from the parent ion. Secondary fragmentation occurs following the release of CO from the initially formed C₃H₇N₂O⁺ channel. This process generates the *m/z* 59.06 (C₂H₇N⁺) fragment, primarily consisting of an ion-molecule complex between neutral methylene imine and the methaniminium cation (m/z 30.03). Additionally, the *m/z* 58.03 (C₂H₄NO⁺) fragment is produced during secondary fragmentation, resulting from the loss of methylene imine and CO neutral species. The methaniminium cation can also be generated via secondary and tertiary fragmentation through the loss of C₂H₃NO + CO or CH₃N + 2 CO from the parent ion. Figure 16a summarizes the time evolution of the different fragmentation products. The theoretical mass spectrum, obtained by summing all the curves leading to the same protonated species, shows relative abundances compared to the

<figure><img src="image_19.png" /><figcaption>FIG. 17. Reaction network of OXA⁰⁰ at a total energy of 178 kcal/mol, detailing the pathways leading to CO loss from the parent ion. The radii of the nodes in this network are proportionally scaled to reflect the self-transition probability of each structure. The blue, gray, and red nodes denote that the nodes belong to the SP, INT, and PF states. The shaded blue circles on the chemical structure highlight the formation of the IMC channel.</figcaption></figure>

primary fragmentation state (see Figure 16b). The spectrum is primarily characterized by the C₃H₇N₂O⁺ species at shorter times and CH₄N⁺ at longer times, consistent with experimental observations where C₃H₇N₂O⁺ dominates at low energy and CH₄N⁺ at higher energies.²³,²⁴

We now summarize the mechanisms leading to the most relevant products obtained from chemical dynamics simulations coupled with graph-theory-based trajectory analysis.

* *Mechanism of Carbon Monoxide Loss, m/z 87.06.* The loss of CO from the parent ion can proceed via two distinct pathways. The first is a direct fragmentation, where the parent ion dissociates into two fragments: a CO molecule and a protonated species. The second pathway involves the formation of an IMC. The mechanisms are identified using the clustering methodology that groups different fragments based on spatial proximity within specific distance thresholds. Figure 16c shows the branching ratio curve for CO loss, distinguishing between direct fragmentation and the IMC pathway. The direct pathway exhibits the highest relative abundance, while the IMC pathway is approximately three times less probable at its peak. The most likely mechanism involves the isomerization of the precursor ion to LIN⁰¹, followed by bond cleavage, resulting in CO release. In Figure 17, we illustrate the fragmentation pathways obtained for OXA⁰⁰ simulations at 178 kcal/mol leading to CO loss. The different structures described in the following are shown in that figure. In particular we detail the fragment-
tation pathways involving OXA<sup>00</sup> (the initial struc-
ture), LIN<sup>01</sup>, and (C<sub>3</sub>H<sub>7</sub>N<sub>2</sub>O<sup>+</sup>)<sub>01</sub><sup>direct</sup> structures. Subse-
quent isomerizations produce various tautomers, such
as (C<sub>3</sub>H<sub>7</sub>N<sub>2</sub>O<sup>+</sup>)<sub>02</sub><sup>direct</sup> and (C<sub>3</sub>H<sub>7</sub>N<sub>2</sub>O<sup>+</sup>)<sub>03</sub><sup>direct</sup>, as well as
cyclic structures like (C<sub>3</sub>H<sub>7</sub>N<sub>2</sub>O<sup>+</sup>)<sub>04</sub><sup>direct</sup>.

The IMC formation mechanism likely begins with the generation of two entities, one being the methaniminium cation, G1. In this process, the neutral component retains its ring structure, followed by a proton transfer to the neutral ring (G2), leading to neutral methylene imine formation. The ring then opens (G3), facilitating CO release via bond cleavage. Ultimately, CO is released, forming a proton-shared complex, (C<sub>3</sub>H<sub>7</sub>N<sub>2</sub>O<sup>+</sup>)<sub>01</sub><sup>imc</sup>, as evidenced by proton transfer dynamics. This mechanism populates several isomers, such as (C<sub>3</sub>H<sub>7</sub>N<sub>2</sub>O<sup>+</sup>)<sub>02</sub><sup>imc</sup> and (C<sub>3</sub>H<sub>7</sub>N<sub>2</sub>O<sup>+</sup>)<sub>03</sub><sup>imc</sup>, which may equilibrate with (C<sub>3</sub>H<sub>7</sub>N<sub>2</sub>O<sup>+</sup>)<sub>01</sub><sup>direct</sup>, indicating an interconnection between direct and IMC pathways. Although this summary addresses the most probable scenarios, the underlying reaction network is more complex.

* *Mechanism of Neutral Methylene Imine Loss, m/z 88.04.* Two mechanisms can be distinguished: direct fragmentation and the IMC pathway, with the direct mechanism being approximately three times more prevalent. The most probable direct pathway involves the release of neutral methylene imine from intermediates (G2) and (G3). In contrast, the IMC pathway involves the formation of (C<sub>3</sub>H<sub>7</sub>N<sub>2</sub>O<sup>+</sup>)<sub>01</sub><sup>imc</sup>. Before CO is released, a proton transfer results in the formyl cation formation, leading to (C<sub>3</sub>H<sub>4</sub>NO<sub>2</sub><sup>+</sup>)<sub>01</sub><sup>imc</sup> and the release of neutral methylene imine. This structure is illustrated in Figure S6b.

* *Mechanism of C<sub>2</sub>H<sub>7</sub>N<sub>2</sub><sup>+</sup> Formation, m/z 59.06.* The formation of C<sub>2</sub>H<sub>7</sub>N<sub>2</sub><sup>+</sup>, involving the release of two CO molecules, appears to occur exclusively via the IMC pathway. This process has a relatively high probability, with a maximum relative abundance of 20% at 5ps. The most probable pathway, as illustrated in Figure S6a, begins with the release of the first CO molecule through the path, OXA<sup>00</sup>, LIN<sup>01</sup>, and (C<sub>3</sub>H<sub>7</sub>N<sub>2</sub>O<sup>+</sup>)<sub>01</sub><sup>direct</sup>, followed by the formation of (C<sub>3</sub>H<sub>7</sub>N<sub>2</sub>O<sup>+</sup>)<sub>03</sub><sup>imc</sup>. From this intermediate, the second CO molecule is released via bond cleavage, leading to the formation of an IMC between the methaniminium cation and neutral methylene imine, ultimately producing (C<sub>2</sub>H<sub>7</sub>N<sub>2</sub><sup>+</sup>)<sub>01</sub><sup>imc</sup>. Notably, several turning points and the formation of different isomers are identified, attributable to proton transfer dynamics.

* *Mechanism of C<sub>2</sub>H<sub>4</sub>NO<sup>+</sup> Formation, m/z 58.03.* The formation of C<sub>2</sub>H<sub>4</sub>NO<sup>+</sup> involves the concurrent loss of a CO molecule and neutral methylene imine. The relative abundance of this channel steadily increases, reaching approximately 10% by the end of the chemical dynamics simulation. Two pathways are possible: direct formation and via the IMC, with both pathways interconnected.

The most probable mechanism follows the same pathway as CO loss up to the formation of (C₃H₇N₂O⁺)ᵢmk. Subsequently, neutral methylene imine is released, leading to C₂H₄NO⁺ formation, as illustrated in Figure S6c (left panel). This species can then dissociate into two components, forming an IMC between neutral methylene imine and a formyl cation, as depicted in Figure S6c (right panel).

* *Mechanism of CH<sub>4</sub>N<sup>+</sup> Formation, m/z 30.03.* The formation of the methaniminium cation is the most abundant channel, becoming predominant after 5ps, see Figure 16b). As previously mentioned, it can be generated via primary, secondary, and tertiary fragmentation, with the primary and tertiary pathways being the most significant. In the case of primary fragmentation, the pathway involves the formation of (G1), where the methaniminium cation is produced. It is important to note that this structure exists in equilibrium with OXA<sup>00</sup> due to the formation of an IMC in the parent channel, IMC<sub>01</sub><sup>01</sup>. The tertiary fragmentation occurs after the dissociation of the IMC created after the loss of two CO molecules. This channel becomes the most abundant after 10ps, see Figure 16a). Finally, the secondary fragmentation pathway involves the breakdown of the IMC created after the loss of one CO molecule, forming (C<sub>3</sub>H<sub>7</sub>N<sub>2</sub>O<sup>+</sup>)<sub>02</sub><sup>imc</sup>.

IV. CONCLUSIONS

Our study reveals significant differences in the fragmenta-
tion dynamics and isomerization pathways between the pro-
tonated oxazolone (OXA) and protonated diketopiperazine
(CYC) systems. This was possible by analyzing an ensem-
ble of reactive trajectories with graph-theory based tools. We
were able thus to compare these detailed analysis on simula-
tions done at the same total energy and with the same amount
of trajectories.

Notably, the number of visited graphs is notably lower in the OXA<sup>00</sup> system (602) compared to the CYC<sup>00</sup> system (1498). Furthermore, the Shannon entropy of the CAN array for OXA<sup>00</sup> (1.77) is lower than that of CYC<sup>00</sup> (2.73), indicating a different entropic tendency toward fragmentation. This difference in entropy may explain why the OXA<sup>00</sup> system exhibits fewer fragmentation channels experimentally, even though it has a lower energy threshold.<sup>24</sup> This observation can be attributed to the dominance of long-lived ion-molecule complexes in OXA<sup>00</sup>, which restricts the exploration of other regions in the potential energy surface before fragmentation occurs. The maximum residence time and clustering methodologies<sup>38</sup> employed in this study have proven to be effective tools for detecting the formation of IMCs.

The OXA<sup>00</sup> system not only presents a lower isomeriza-
tion and reaction energy threshold than its CYC<sup>00</sup> counter-
part but also exhibits a reaction threshold difference of 0.9 eV,
consistent with experimental findings.<sup>24</sup> Although our simula-
tions cannot be directly compared to experimental results due
to the differences in time scales by several orders of magni-
tude, they provide valuable insight when comparing the behavior of systems under the same conditions. The isomerization landscape of OXA⁰⁰ is largely governed by the formation of IMCs, in contrast to the CYC⁰⁰ system, where linear and cyclic structures are more prominent. Furthermore, the isomerization pathway between them implicates the formation of the IMCs intermediates as a critical step. Our graph-theory-based methodology effectively elucidated the formation pathways of these IMCs, demonstrating the utility of identifying turning points as a strategy for detecting IMCs.

## AUTHOR CONTRIBUTIONS

A.F.P.M., T.B., and R.S. conceptualized the project. A.F.P.M. performed the simulations and data analysis. T.B. and R.S. provided the computational resources. A.F.P.M. prepared the first draft of the manuscript with input from all the authors.

## CONFLICTS OF INTEREST

The authors declare no competing interests.

## ACKNOWLEDGEMENTS

This work was partially supported by ANR DynBioReact (Grant No. ANR-14-CE06-0029-01). A.F.P.M. and T.B. thank the University of Geneva for the financial support. We thank MsC. Elena Zerbato, Dr. Arnulf Rosspeintner, Dr. Max Latevi, and Dr. Alejandro Gutierrez Quintanilla for useful discussions and suggestions on the manuscript.

## DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## BIBLIOGRAPHY

<footnote n="1">¹I. A. Papayannopoulos, "The interpretation of collision-induced dissociation tandem mass spectra of peptides," Mass Spectrom. Rev. **14**, 49–73 (1995).</footnote>

<footnote n="2">²J. Zaia, "Mass spectrometry of oligosaccharides," Mass Spectrom. Rev. **23**, 161–227 (2004).</footnote>

<footnote n="3">³B. Paizs and S. Suhai, "Fragmentation pathways of protonated peptides," Mass Spectrom. Rev. **24**, 508–548 (2004).</footnote>

<footnote n="4">⁴J. Laskin and J. H. Futrell, "Collisional activation of peptide ions in ft-icr mass spectrometry," Mass Spectrom. Rev. **22**, 158–181 (2003).</footnote>

<footnote n="5">⁵X. Han, A. Aslanian, and J. R. Yates III, "Mass spectrometry for proteomics," Curr. Opin. Chem. Biol. **12**, 483–490 (2008).</footnote>

<footnote n="6">⁶F. W. McLafferty, K. Breuker, M. Ji, X. Han, G. Infusini, H. Jiang, X. Kong, and T. P. Begley, "Top-down ms, a powerful complement to the high capabilities of proteolysis proteomics," FEBS J. **274**, 6256–6268 (2007).</footnote>

<footnote n="7">⁷P. Roepstorff and J. Fohlman, "Proposal for a common nomenclature for sequence ions in mass spectra of peptides," Biomed. Mass Spectrom. **11**, 601 (1984).</footnote>

<footnote n="8">⁸K. Biemann, "Contributions of mass spectrometry to peptide and protein structure," Biomedical & Environmental Mass Spectrometry **16**, 99–111 (1988), https://onlinelibrary.wiley.com/doi/pdf/10.1002/bms.1200160119.</footnote>

<footnote n="9">⁹A. G. Harrison, A. Young, C. Bleiholder, S. Suhai, and B. Paizs, "Scrambling of sequence information in collision-induced dissociation of peptides," J. Am. Chem. Soc. **128**, 10364–10365 (2006).</footnote>

<footnote n="10">¹⁰A. G. Harrison, "To b or not to b: The ongoing saga of peptide b ions," Mass Spectrom. Rev. **28**, 640–654 (2009).</footnote>

<footnote n="11">¹¹I. S. Saminathan, X. S. Wang, Y. Guo, O. Krakovska, S. Voisin, A. C. Hopkinson, and K. M. Siu, "The extent and effects of peptide sequence scrambling via formation of macrocyclic b ions in model proteins," J. Am. Soc. Mass Spectrom. **21**, 2085–2094 (2010).</footnote>

<footnote n="12">¹²M. Tirado and N. C. Polfer, "Defying entropy: forming large head-to-tail macrocycles in the gas phase," Angew. Chem. Int. Ed. **51**, 6436–6438 (2012).</footnote>

<footnote n="13">¹³T. Yalcin, C. Khouw, I. G. Csizmadia, M. R. Peterson, and A. G. Harrison, "Why are b ions stable species in peptide spectra?" J. Am. Soc. Mass Spectrom. **6**, 1165–1174 (1995), pMID: 24214067, https://doi.org/10.1016/S1044-0305(95)00569-2.</footnote>

<footnote n="14">¹⁴K. Eckart, M. C. Holthausen, W. Koch, and J. Spiess, "Mass spectrometric and quantum mechanical analysis of gas-phase formation, structure, and decomposition of various b2 ions and their specifically deuterated analogs," J. Am. Soc. Mass Spectrom. **9**, 1002–1011 (1998), https://doi.org/10.1016/S1044-0305(98)00076-2.</footnote>

<footnote n="15">¹⁵A. G. Harrison, I. G. Csizmadia, and T.-H. Tang, "Structure and fragmentation of b2 ions in peptide mass spectra," J. Am. Soc. Mass Spectrom. **11**, 427–436 (2000), pMID: 10790847, https://doi.org/10.1016/S1044-0305(00)00104-5.</footnote>

<footnote n="16">¹⁶B. Paizs, S. Suhai, and A. G. Harrison, "Experimental and theoretical investigation of the main fragmentation pathways of protonated h-gly-gly-sar-oh and h-gly-sar-sar-oh," J. Am. Soc. Mass Spectrom. **14**, 1454–1469 (2003), pMID: 14652192, https://doi.org/10.1016/j.jasms.2003.07.001.</footnote>

<footnote n="17">¹⁷N. C. Polfer and J. Oomens, "Reaction products in mass spectrometry elucidated with infrared spectroscopy," Phys. Chem. Chem. Phys. **9**, 3804–3817 (2007).</footnote>

<footnote n="18">¹⁸R. K. Sinha, U. Erlekam, B. J. Bythell, B. Paizs, and P. Maître, "Diagnosing the protonation site of b2 peptide fragment ions using irmpd in the x–h (x = o, n, and c) stretching region," J. Am. Soc. Mass Spectrom. **22**, 1645–1650 (2011).</footnote>

<footnote n="19">¹⁹B. R. Perkins, J. Chamot-Rooke, S. H. Yoon, A. C. Gucinski, Somogyi, and V. H. Wysocki, "Evidence of diketopiperazine and oxazolone structures for ha b2+ ion," J. Am. Chem. Soc. **131**, 17528–17529 (2009), pMID: 19904982, https://doi.org/10.1021/ja9054542.</footnote>

<footnote n="20">²⁰S. Zou, J. Oomens, and N. C. Polfer, "Competition between diketopiperazine and oxazolone formation in water loss products from protonated arggly and glyarg," Int. J. Mass Spectrom. **316-318**, 12–17 (2012), special Issue Honoring Alex G. Harrison on the Occasion of his 80th Birthday: Structures, Reactions, and Thermochemistry of Gas-Phase Ions.</footnote>

<footnote n="21">²¹J. S. Klassen and P. Kebarle, "Collision-induced dissociation threshold energies of protonated glycine, glycinamide, and some related small peptides and peptide amino amides," J. Am. Chem. Soc. **119**, 6552–6563 (1997), https://doi.org/10.1021/ja962813m.</footnote>

<footnote n="22">²²C. F. Rodriquez, A. Cunje, T. Shoeib, I. K. Chu, A. C. Hopkinson, and K. W. M. Siu, "Proton migration and tautomerism in protonated triglycine," J. Am. Chem. Soc. **123**, 3006–3012 (2001), pMID: 11457011, https://doi.org/10.1021/ja0015904.</footnote>

<footnote n="23">²³H. E. Aribi, C. F. Rodriquez, D. R. P. Almeida, Y. Ling, W. W.-N. Mak, A. C. Hopkinson, and K. W. M. Siu, "Elucidation of fragmentation mechanisms of protonated peptide ions and their products: a case study on glycylglycylglycine using density functional theory and threshold collision-induced dissociation," J. Am. Chem. Soc. **125**, 9229–9236 (2003), pMID: 15369379, https://doi.org/10.1021/ja0207293.</footnote>

<footnote n="24">²⁴P. B. Armentrout and A. A. Clark, "The simples $b_2^+$ ion: Determining its structure from its energetics by a direct comparison of the threshold collision-induced dissociation of protonated oxazolone and diketopiperazine," Phys. Chem. Chem. Phys. **316-318**, 182–191 (2012).</footnote>
25 W. L. Hase, "Simulations of gas-phase chemical reactions: Applications to s<sub>sub</sub>n<sub>n</sub>/<sub>sub</sub>2 nucleophilic substitution," *Science* **266**, 998–1002 (1994), https://www.science.org/doi/pdf/10.1126/science.266.5187.998.

26 Z. Homayoon, V. Macaluso, A. Martin-Somer, M. C. N. B. Muniz, I. Borges, W. L. Hase, and R. Spezia, "Chemical dynamics simulations of cid of peptide ions: comparisons between tik(h+)2 and tlk(h+)2 fragmentation dynamics, and with thermal simulations," *Phys. Chem. Chem. Phys.* **20**, 3614–3629 (2018).

27 A. Pérez-Mellor, I. Alata, V. Lepere, R. Spezia, and A. Zehnacker-Rentien, "Stereospecific collision-induced dissociation and vibrational spectroscopy of protonated cyclo (tyr-pro)," *Int. J. Mass Spectrom.* **465**, 116590 (2021).

28 Pérez-Mellor, Ariel, Le Barbu-Debus, Katia, Lepere, Valeria, Alata, Ivan, Spezia, Riccardo, and Zehnacker, Anne, "Structure and collision-induced dissociation of the protonated cyclo his-phe dipeptide: mechanistic studies and stereochemical effects," *Eur. Phys. J. D* **75**, 165 (2021).

29 G. H. Peslherbe and W. L. Hase, "Semiempirical MNDO, AM1, and PM3 direct dynamics trajectory studies of formaldehyde unimolecular dissociation," *J. Chem. Phys.* **104**, 7882–7894 (1996), https://pubs.aip.org/aip/jcp/article-pdf/104/20/7882/19011839/7882_1_online.pdf.

30 A. Malik, R. Spezia, and W. L. Hase, "Unimolecular fragmentation properties of thermometer ions from chemical dynamics simulations," *J. Am. Soc. Mass Spectrom.* **32**, 169–179 (2021), pMID: 33210535, https://doi.org/10.1021/jasms.0c00200.

31 I. A. Hubner, E. J. Deeds, and E. I. Shakhnovich, "Understanding ensemble protein folding at atomic detail," *Proc. Natl. Acad. Sci. (USA)* **103**, 17747–17752 (2006), https://www.pnas.org/doi/pdf/10.1073/pnas.0605580103.

32 K.-C. Chou, "Applications of graph theory to enzyme kinetics and protein folding kinetics: Steady and non-steady-state systems," *Biophys. Chem.* **35**, 1–24 (1990).

33 J. Koča, S. Pérez, and A. Imberty, "Conformational analysis and flexibility of carbohydrates using the cicada approach with mm3," *J. Comput. Chem.* **16**, 296–310 (1995), https://onlinelibrary.wiley.com/doi/pdf/10.1002/jcc.540160305.

34 S. Bougueroua, R. Spezia, S. Pezzotti, S. Vial, F. Quessette, D. Barth, and M.-P. Gaigeot, "Graph theory for automatic structural recognition in molecular dynamics simulations," *J. Chem. Phys.* **149**, 184102 (2018), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.5045818/13232918/184102_1_online.pdf.

35 O. M. Becker and M. Karplus, "The topology of multidimensional potential energy surfaces: Theory and application to peptide structure and kinetics," *J. Chem. Phys.* **106**, 1495–1517 (1997), https://pubs.aip.org/aip/jcp/article-pdf/106/4/1495/19148343/1495_1_online.pdf.

36 E. Martínez-Núñez, "An automated transition state search using classical trajectories initialized at multiple minima," *Phys. Chem. Chem. Phys.* **17**, 14912–14921 (2015).

37 D. R. Galimberti, S. Bougueroua, J. Mahé, M. Tommasini, A. M. Rijs, and M.-P. Gaigeot, "Conformational assignment of gas phase peptides and their h-bonded complexes using far-ir/thz: Ir-uv ion dip experiment, dft-md spectroscopy, and graph theory for mode assignment," *Faraday Discuss.* **217**, 67–97 (2019).

38 A. F. Perez-Mellor and R. Spezia, "Determination of kinetic properties in unimolecular dissociation of complex systems from graph theory based analysis of an ensemble of reactive trajectories," *J. Chem. Phys.* **155**, 124103 (2021), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/5.0058382/15616903/124103_1_online.pdf.

39 D. Townsend, S. A. Lahankar, S. K. Lee, S. D. Chambreau, A. G. Suits, X. Zhang, J. Rheinecker, L. B. Harding, and J. M. Bowman, *Science* **306**, 1158–1161 (2004).

40 J. M. Bowman and P. L. Houston, "Theories and simulations of roaming," *Chem. Soc. Rev.* **46**, 7615–7624 (2017).

41 V. Macaluso, D. Scuderi, M. E. Crestoni, S. Fornarini, D. Corinti, E. Dalloz, E. Martinez-Nunez, W. L. Hase, and R. Spezia, "l-cysteine modified by s-sulfation: Consequence on fragmentation processes elucidated by tandem mass spectrometry and chemical dynamics simulations," *J. Phys. Chem. A* **123**, 3685–3696 (2019), pMID: 30945861, https://doi.org/10.1021/acs.jpca.9b01779.

42 G. B. Rocha, R. O. Freire, A. M. Simas, and J. J. P. Stewart, "Rm1: A reparameterization of am1 for h, c, n, o, p, s, f, cl, br, and i," *J. Comput. Chem.* **27**, 1101–1111 (2006), https://onlinelibrary.wiley.com/doi/pdf/10.1002/jcc.20425.

43 S. Grimme, "Accurate description of van der waals complexes by density functional theory including empirical corrections," *J. Comput. Chem.* **25**, 1463–1473 (2004), https://onlinelibrary.wiley.com/doi/pdf/10.1002/jcc.20078.

44 G. L. Barnes and W. L. Hase, "Energy transfer, unfolding, and fragmentation dynamics in collisions of n-protonated octaglycine with an h-sam surface," *J. Am. Chem. Soc.* **131**, 17185–17193 (2009), pMID: 19929018, https://doi.org/10.1021/ja904925p.

45 Z. Homayoon, S. Pratihar, E. Dratz, R. Snider, R. Spezia, G. L. Barnes, V. Macaluso, A. Martin Somer, and W. L. Hase, "Model simulations of the thermal dissociation of the tik(h+)2 tripeptide: Mechanisms and kinetic parameters," *J. Phys. Chem. A* **120**, 8211–8227 (2016).

46 G. L. Barnes, K. Young, L. Yang, and W. L. Hase, "Fragmentation and reactivity in collisions of protonated diglycine with chemically modified perfluorinated alkylthiolate-self-assembled monolayer surfaces," *J. Chem. Phys.* **134**, 094106 (2011), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.3558736/13306826/094106_1_online.pdf.

47 W. L. Hase and D. G. Buckowski, "Monte carlo sampling of a microcanonical ensemble of classical harmonic oscillators," *Chem. Phys. Lett.* **74**, 284–287 (1980).

48 W. L. Hase, R. J. Duchovic, X. Hu, A. Komornicki, K. F. Lim, D.-H. Lu, G. H. Peslherbe, K. N. Swamy, S. R. V. Linde, A. Varandas, H. Wang, and R. J. Wolf, "Venus, a general chemical dynamics computer program," *QCPE Bull.* **16**, 671 (1996).

49 J. J. P. Stewart, L. K. Fiedler, J. Zheng, I. Rossi, W.-P. Hu, G. C. Lynch, Y.-P. Liu, P. Zhang, Y.-Y. Chuang, and J. P. et al., "MOPAC, version 5.022mn based on MOPAC 5.0," (2015).

50 B. D. McKay and A. Piperno, "Practical graph isomorphism, ii," *J. Symbolic Comput.* **60**, 94 – 112 (2014).

51 S. Even, *Graph Algorithms*, 2nd ed., edited by G. Even (Cambridge University Press, 2011).

52 E. W. Dijkstra, "A note on two problems in connexion with graphs," *Numerische Mathematik* **1**, 269–271 (1959).

53 R. F. Grote and J. T. Hynes, "," *J. Chem. Phys.* **73**, 2715 (1980).

54 T. Baer and W. L. Hase, *Unimolecular reaction dynamics: theory and experiments*, (Oxford University Press, New York, 1996).

55 I. N. Levine, *Quantum Chemistry*, 7th ed. (Pearson, 2014).

56 M. Brehm and B. Kirchner, "Travis - a free analyzer and visualizer for monte carlo and molecular dynamics trajectories," *J. Chem. Info. Model.* **51**, 2007–2023 (2011), pMID: 21761915, https://doi.org/10.1021/ci200217w.

57 N. M. Rodiguin and E. H. Rodiguina, *Consecutive Chemical Reactions: Mathematical Analysis and Development*, (D. Van Nostrand Co., New York, 1964).

58 R. Spezia and H. Dammak, "On the use of quantum thermal bath in unimolecular fragmentation simulation," *J. Phys. Chem. A* **123**, 8542–8551 (2019), pMID: 31509415, https://doi.org/10.1021/acs.jpca.9b06795.