Ariel F Perez-Mellor

Riccardo Spezia

Determination of kinetic properties in unimolecular dissociation of complex systems from graph-theory based analysis of an ensemble of reactive trajectories

Keywords:

teaching and research institutions in France or abroad, or from public or private research centers.

Determination of kinetic properties in unimolecular dissociation of complex systems from graph-theory based analysis of an ensemble of reactive trajectories

I. INTRODUCTION

Using dynamical simulations to understand chemical reactivity represents a wide and important field in physical chemistry. [1][2][3][4][5] More recently, thanks to the improved computer power and new developments in electronic structure theory, it becomes possible to follow on-the-fly bond breaking/making processes also for relatively complex molecules, like peptides, sugars or polycyclic aromatic hydrocarbons. [6][7][8] While fitted potential energy surfaces are often used to understand in detail reactions of relatively small systems, [9][10][11][12][13] the coupling between dynamics and electronic structure theory made it possible to study the evolution on-the-fly of relatively large systems. For example, uni-and bi-molecular reaction dynamics simulations were largely studied with the aim of understanding several processes, from combustion to atmospheric chemistry, from interfacial reactions to photolysis or astrochemistry. [14][15][16][17][18][19][20] Unimolecular reactivity represents an important class of reactions, [START_REF] Baer | Unimolecular reaction dynamics: theory and experiments[END_REF] with a relevant application in the field of mass spectrometry. [START_REF] Song | Theoretical Mass Spectrometry[END_REF] Recently, trajectories-based methods were used to understand several gas-phase ion chemistry experiments, like electron ionization mass spectrometry, [START_REF] Grimme | [END_REF][24][25] surface-induced dissociation (SID) 6,26,27 or collision-induced dissociation. 28,29 a) ariel.perezmellor@unige.ch; Present address: Department of Physical Chemistry, University of Geneva, 30 Quai Ernest-Ansermet, 1211 Geneva 4, Switzerland. b) riccardo.spezia@sorbonne-universite.fr

In particular, in the field of tandem mass spectrometry, simulations were able to predict fragmentation products, propose reaction mechanisms and quantify reaction energy thresholds for both statistical and non-statistical processes. 7 Statistical unimolecular fragmentation kinetics of relatively large systems was obtained by activating the initial system with excess vibrational energy. By following the initial product's decay as a function of time, it was possible to obtain unimolecular rate constants and, from Arrhenius-like plots, the corresponding activation energies. 16,[30][31][32] On-the-fly simulations are based on the time propagation of atomic positions expressed in cartesian coordinates. Thus, from the computational point of view, a molecule is a set of atoms connected by chemical bonds. However, since the system is represented by interacting nuclei and electrons, the bonds are not directly defined. One can state the existence and evolution of such bonds by geometrical and/or electron density criteria. When systems are simple and trajectories relatively few, the outcome of simulations can be analyzed either by ad hoc quantifications (distances, angles, etc ...) or by visual inspection of the resulting trajectories. However, thanks to computational improvements over the last years, it is now possible to run several and long trajectories of relatively big systems for which different and complex behaviors are possible. Of course, it is impossible to predict in advance all the possible reaction pathways only based on chemical intuition. The results, however, consist of a huge amount of data from which advanced analysis tools, like graph-theory, are needed to obtain physical properties. For this reason, graph-theory-based methods were used and developed to better understand the behavior of complex molecular systems from a relatively important amount of data. Examples can be found in protein folding [START_REF] Hubner | Proc. Natl. Acad. USA[END_REF] , enzyme kinetics [START_REF] Chou | [END_REF] , conformational analysis 35 , protein structure and flexibility identification 36,37 , peptide structure and kinetics studies, 38 or assignment of vibrational normal modes. 39 Molecular dynamics simulations are widely used in different fields, and they provide a huge amount of data from which important molecular properties should be obtained. To this end, several authors used graph theory to determine the different structures sampled during long simulations. [40][41][42][43][44] Pietrucci and Andreoni developed another useful application to drive biased molecular dynamics simulations from reactants to products without imposing elaborated geometrical constraints. 45,46 Recently, an original way to reduce computational time when analyzing a large set of data issued from molecular simulations using graph theory was proposed by Bougueroua et al. 47 Martinez-Nunez and co-workers proposed a specific application to chemical reactivity. [48][49][50] It is based on molecular dynamics in which the graph theory algorithms automatically determine minima and saddle points along a reaction path. This method was applied to photo-fragmentation, 51 unimolecular fragmentation, 52 and organometallic catalysis. 53 We recently used graph theory to automatically detect if a molecular system dissociates after internal energy activation and obtain, from an ensemble of quasi-classical trajectories, the distribution of molecular fragments. 54 From this information, an in silico mass spectrum can be obtained as well as unimolecular rate constants. 30,55 However, different steps were determined ad hoc or not fully automatized.

In the present work, we show how graph-theory can be developed to obtain the structural and kinetic information needed to describe the reactivity of complex molecular systems as obtained from an ensemble of trajectories, here focusing in particular on gas phase unimolecular fragmentation. In particular, we consider the possibility that an initial structure passes through intermediates before dissociating. One important point is to correctly consider the possibility that the system oscillates between intermediate(s) and reactant before it eventually dissociates. In mass spectrometry, it is well known that ion-molecule complexes can be formed before fragmentation. 56,57 A recent study of L-cysteine fragmentation 52 has shown that this can correspond to the dynamical mechanisms called roaming. 58,59 The accurate detection of such ion-molecule complexes' formation can be problematic due to their elusive dynamical nature. The specific analysis of trajectories to identify such states is discussed here.

The methodology is explained in general forms and applied to the prototypical case of fragmentation of protonated cyclo-di-glycine, for which experimental mass spectra are available and fragments known. 60 It will serve as an example to show how we can obtain from graph-theory structural and kinetic information of fragmentation of a large system. In particular, peptides fragmentation represents a typical example in which the correct identification of intermediate structures is crucial to explain reactivity: after activation, the excess proton often moves to different basic sites, weakening the adjacent bonds and inducing fragmentation. This is the so-called "mobile proton model," which was shown to be a general mechanism typical of peptide fragmentation [61][62][63] but also relevant in other compounds, like nucleic acids. 55,64 In recent works, we have successfully applied some aspects of the methodology proposed here to qualitatively disentangle the complex dynamics in CID experiments of dipeptides that contain the DKP motif, such as cyclo Tyr-ProH +65 and cyclo Phe-HisH + . 66 This work illustrates how graph theory can be applied to quantitatively describe such complex unimolecular processes, where the protonated cyclo-di-glycine is used as a test case. In setion II, we first explain how chemical dynamics simulations are performed: they provide this specific example's trajectories, but the methodology is more general and does not depend on how simulations are performed. In section III, we show how graph theory formalism can be adapted to analyze an ensemble of trajectories. We have then applied it to a simple three-state kinetic model (section IV) from which overall properties can be obtained. In section V, we show how the model can be extended to analyze in more details one specific pathway. Finally, we demonstrate how the method can be used to obtain a mass spectrum. Section VII concludes the article.

II. CHEMICAL DYNAMICS SIMULATIONS

Two isomers of cyclo-di-glycine were used as initial structures for the dynamical simulations: the most stable cyclic structure (labeled CYC 00 ) and a linear isomer (LIN 00 ) which is higher in energy and largely populated. They are shown in Figure 1.

CYC 00 LIN 00 Trajectories were propagated on-the-fly using the semiempirical Hamiltonian RM1-D 67,68 to obtain energies and gradients of the two systems. Note that this method was shown to correctly characterize dissociation properties of protonated peptides in the gas phase 31,69,70 and includes dispersion corrections.

To induce unimolecular fragmentations, the systems were activated via excess vibrational energy through microcanonical normal mode sampling. 71 We used several internal energy values in the 131-217 kcal/mol range, listed in Table I (labeled as E v ). Rotational energy was incorporated via a Boltzmann distribution at 300 K. For each system and total energy, an ensemble of trajectories (about 10 000 per point) was run, integrating Newton's equations of motion numerically via the velocity Verlet algorithm with a time step of 0.1 fs, ensuring energy conservation. Each trajectory was propagated up to 20 or 40 ps depending on the energy. Details on simulation conditions are summarized in Table I, where Λ is the simulation time of each trajectory. In the same table we also report the average potential ( V ) and kinetic ( K ) energies (with corresponding standard deviations) of the activated structures obtained from simulations with the energy decomposition method discussed in section IV B where they will also be discussed. The zero of the potential energy was set at the electronic energy of CYC 00 . Table I. Simulation details. Λ is the total integration time, N the number of trajectories, E T , the total energy, E v the initial excess vibrational energy, K the average kinetic energy and V the average potential energy, where σ K and σ V are the corresponding standard deviations.

Λ N E T E v K σ K V σ V ps kcal/
All simulations were performed with the general chemical dynamics software Venus 72 which is interfaced with MOPAC 5.022mn [START_REF] Stewart | MOPAC, version 5.022mn based on MOPAC 5[END_REF] for electronic structure calculations.

III. GRAPH THEORY ANALYSIS OF AN ENSEMBLE OF REACTIVE TRAJECTORIES

Chemical dynamics trajectories provide a huge amount of data in terms of positions (and momenta) of the different atoms as a function of time. We now define the key quantities we can obtain by applying graph theory to this set of data. These quantities will be used to obtain the different physical quantities. In the present application we will focus in particular on what is relevant to obtain rate constants and (time dependent) theoretical mass spectrum.

A. Data storage

In all trajectories, the information is saved every τ steps. In the present case we have chosen τ = 50 fs, and this will depend on the characteristic time-scale of the process of interest and on the storage capabilities. As a result of this time discretization, several 2D-arrays are constructed to store and analyze the information appropriately.

In the following, X [i, j] denotes a generic 2D-array where i indicates the trajectory number (i ∈ [1, N]) and j the snapshot (j ∈ [1, M] with M = Λ/τ + 1). We have, for example, the atomic positions XYZ [i, j], the potential POT [i, j], the kinetic KIN [i, j], and total energy of the system TOE [i, j].

In general, two kinds of arrays X [i, j] can be distinguished according to the type of variable it might contain: continuous or discrete. X [i, j] is said to be a continuous (discrete) 2D-array if and only if (iff) the elements X ij are continuous (discrete). A continuous array can be transformed into a discrete array through a discretization process, which will consist of introducing threshold values. A typical example used in this work is transforming the cartesian XYZ [i, j] array into the adjacency ADJ [i, j] array, as explained in the next subsection.

B. Building undirected and unweighted κ-Coloured Graph

A critical aspect of the unimolecular fragmentation simulation is to identify whether an isomerization or fragmentation event occurs from a computer analysis of time evolution of cartesian coordinates without any visual inspection. At this end we used different tools of graph theory. [START_REF] Trinajstic | Chemical Graph Theory[END_REF] The analysis of a set of N-trajectories begins with the treatment of the aforementioned XYZ [i, j]-array. It is a 2D-array where each XYZ ij element contains the position of the atoms for each trajectory i at the discretized time (j -1) τ.

The geometry contained in XYZ ij can be converted into a simple undirected and unweighted κ-colored graph G = (V, E) by following a discretization criterion that establishes whether two atoms are linked or not. The finite set of vertices V = {v 1 , . . . , v λ } denotes all the λ -atoms while the finite set E ⊆ V × V of edges (v n , v m ) with m, n ∈ [1, λ ] indicates all the possible connections between them. Each vertex v n is colored with a suitable color C = {1, 2, . . . κ} according to the type of atom it represents. Therefore, the maximum number of different types of chemical atoms in the molecule defines the color number κ of G. In addition, the coloring of the vertices π ∈ V × C is known as a partition of the nodes. These concepts are directly taken from standard graph theory. [START_REF] Trinajstic | Chemical Graph Theory[END_REF] A practical way of representing a graph G is through its adjacency matrix A ∈ R λ ×λ where λ = |V|, the number of vertices corresponding to the number of atoms. It is a symmetric (0, 1)-matrix in which the diagonal elements are all zeros. The rest of the elements A nm are one iff exists an edge between vertex v n and v m , otherwise they are zero. To build the A matrix for a given XYZ ij molecular geometry we define the threshold bond distances as the sum of the atomic radii of the two atoms R n + R m involved in the bond, multi-plied by a factor α:

A nm =    1 →| r n -r m |≤ α • (R n + R m ) 0 →| r n -r m |> α • (R n + R m )
The α parameter allows for bond elongation during the molecular dynamics, i.e., vibrations around stables geometries. It means that all the geometries around stable isomers have the same G (V, E). The optimum α value of 1.264 was determined in previous studies 30,54 since it is able to provide reasonable chemical structures avoiding unphysical atomic valences. For the atomic radii, we used previously reported values 54 , and notably: The α value remains constant in all the systems studied in this work. To overcome the dependency of α with vibrational energy, some corrections were introduced during the definition of the states, as will be detailed in the section IV A.

Thus, for each XYZ ij , there exists a simple undirected κ-coloured graph that can be represented by its adjacency matrix ADJ ij . The partition is built with the set of atom types respecting the order of appearance in the XYZ file; in our case a possible representation is thus π = {HHHHHHHCCCCNNOO}. This representation was chosen to address events such as isomerization and reactivity specifically. Simultaneously, structures that differ, for example, by the rotation of one or more atoms or groups, will not be considered different. Shape descriptors such as inertia moments or radius of gyration could be used to discriminate between rotamers, for example, but it is out of the scope in this work which is devoted on describing chemical reactivity and thus not considered presently.

In Figure 2, we show a snapshot extracted from the simulations to illustrate all the steps of the methodology developed here. Table II displays the adjacency matrix of the structure sketched in Figure 2 constructed following the steps mentioned above.

C. Derived Information From ADJ [i, j]-array

The analysis of the 2D-array ADJ [i, j] through graphtheory tools is a key point of this study. In the following we detail some useful 2D-arrays that can be obtained from

ADJ [i, j]. • ADJ [i, j] + π → CAN [i, j].
Having defined the partition π of the graph ADJ ij and by using a graph isomorphism algorithm, one can get the canonical label, CAN ij , of ADJ ij . Therefore, each CAN ij -element is the canonical label of its corresponding ADJ ij . The labelling process was carried out by the use of the amtog and labelg tools included in the NAUTY package [START_REF] Mckay | [END_REF] . The introduction of this graph invariant allows to overcome the isomorphism problem which otherwise will occur in the automatic analysis of the trajectories. H1 H2 H3 H4 H5 H6 H7 C1 C2 C3 C4 N1 N2 O1 O2 H1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 H2 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 H3 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 H4 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 H5 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 H6 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 H7 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 C1 1 0 1 0 0 0 0 0 0 0 0 1 0 0 0 C2 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 C3 0 0 0 0 0 0 1 0 0 0 0 0 1 1 0 C4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 N1 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 N2 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 O1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 O2 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 Table II. An adjacency matrix of the snapshot depicted in Figure 2. The bonds were colored accordingly to the fragment to which they belong.

Therefore, one can count precisely the amount of different graphs that are formed during the dynamics. This precise counting is impossible using only the adjacency matrix formalism. The permutation of two rows or two columns of the same atom type inherently implies a change in the adjacency matrix. At the same time, the resulting structure is chemically identical to that of the precursor. A typical example is when two hydrogen atoms exchange themselves. In other terms, canonical labeling is an efficient way to identify if two graphs are identical or not. Two graphs are said to be identical iff they share the same canonical label; otherwise, they are different. On the top of Figure 2 is represented the canonical label (CL) of the adjacency matrix of the entire snapshot (Table II) by using the partition introduced above. The Graph6 format was used for the canonical representation. [START_REF] Mckay | [END_REF] • ADJ [i, j] → CLU [i, j]. Fruitful information can be extracted by taking into account the connectivity of the graph. The graph ADJ ij can be classified as connected or disconnected. A graph is said to be connected when every pair of vertices in the graph is connected, it means that there is a path between v n and v m . Chemically speaking, this corresponds to an undissociated molecule. Otherwise, they are called disconnected, which corresponds to two or more molecules. Moreover, every disconnected graph is made by a set of connected subgraphs. By identifying and counting all the connected subgraphs in ADJ ij , using the depth-first search (DFS) algorithm [START_REF] Even | Graph Algorithms[END_REF] , it is possible to decompose ADJ ij into all the connected entities, such that,

ADJ ij = D ∑ k=1 ADJ k ij (1)
where, D is the number of disconnected components of the associated graph ADJ ij . An example of this decomposition is given in Table II, where the full matrix can be seen as the sum of the matrices of the three resulting fragments, each labeled with a different color (red, blue and green). Then, the canonical label of every ADJ k ij is calculated by using the same partition π previously defined (labeled CL1, CL2 and CL3 in Figure 2). The latter allows an adequate count of the different fragments.

Having determined the atoms that belong to a given fragment, we calculate (in the case of a disconnected graph) the centers of mass (COM) of every connected subgraph and the distances between them (labeled d12, d23, and d13 in Figure 2). This is an important information for determining the formation of ionmolecule complexes as explained further below.

All this information is stored in the CLU [i, j] array.

• ADJ [i, j] → CLU [i, j] + M → MAS [i, j].
• ADJ [i, j] → CLU [i, j] + M → MAS [i, j] → CMAS [i, j].
Similar to the α parameter, we introduce a cut-off distance (dclus) between the COM to determine whether two connected subgraphs form a complex or not. Then, by inspecting the corresponding adjacency matrix B ∈ R D×D , we can detect such events. The B matrix is defined as,

B pq =    1 →| R CM p -R CM q |≤ dclus 0 →| R CM p -R CM q |> dclus
where p and q run from one up to the number of disconnected graphs identified in the snapshot. Once the formation of an ion-complex is detected, their masses are calculated and stored in CMAS. As a consequence, each CMAS ij -element corresponds to the set of cluster mass recognized in ADJ ij . For example, CMAS ij = {M1 + M2, M3} if a complex is formed between M1 and M2 as in the example of Figure 2.

This approach of using only the distance dclus as a criterion to determine whether two fragments form a complex or not could be not suitable for larger systems assuming extended configurations (as can occur in peptides and more in general in large polymers). This is due to the fact that the distance between the outer atom in the fragment and the COM can be considerably large. A possible solution can be to divide the whole molecule in sub-units and then perform the same analysi with respect to the COM of each subunit. In the case of (bio)-polymers, one can use as sub-units the different monomers which are normally relatively small and compact.

The dclus parameter was obtained from the analysis of trajectories which form disconnected graphs, and in particular the CYC 00 simulations at 177 kcal/mol and LIN 00 at 151 kcal/mol. For the two sets we calculate the distances between the COM of the different disconnected graphs obtained from which we get a pair distribution function of fragment distances. The normalized distribution functions (PDF) are shown in Figure 3. The two PDFs show a maximum at around 3.75 Å which corresponds mainly to the transient separation of hydrogen atoms from the main structure as detected by the algorithm (this does not correspond to a physical phenomenon of dissociation, but to a hot vibrational state), and sketched in Figure 3. Thus, this 3.57 Å distance can be seen as an effective radius of all the connected graphs, and it can serve as a reference to establish the cut-off for ion-molecule complex detection, which must be larger than it. A cut-off of 8 Å is more than twice this value and, as shown in Figure 3 can represent a good estimation of a distance to automatically identify ion-molecule complexes from the analysis of an ensemble of trajectories. of the associated graph, respectively. The number of edges can be computed directly from the adjacency matrix, such that:

• ADJ [i, j] → CLU [i, j] → CYC [i, j].
|E| = 1 2 |V| ∑ n=1 |V| ∑ m=1 A nm (2)
while the value of D was already obtained in CLU [i, j]. Therefore, the CYC ij -element account for the number of cycles contained in the adjacency matrix ADJ ij . For example, if c = 1 the graph (and so the corresponding molecule) has one ring.

Table III summarizes all the information described and extracted from the simulations. With all of them in our hands, it is possible to identify several events occurring in the simulations, such as isomerization, fragmentation, or ion-molecule complex formation.

Before entering into the detailed analysis of the trajectories, we will present some properties of the 2D arrays that will be essential for the analysis and further discussion.

D. Properties of a discrete X [i, j]
All the discrete X [i, j] 2D-arrays described previously fulfill the following properties:

1. Maximum number of different items. X [i, j] has a maximum number of different items ρ X ∈ [1, N • M].
The minimum value of ρ X corresponds to one if all the elements of X [i, j] are identical while its maximum value corresponds to N • M if all the items are different. Therefore, the sample space of X can be defined Table III. Summary of all the information derived and stored from the present analysis conducted on the ensemble of trajectories.

as Ω X ∈ {X 1 , X 2 , ..., X ρ }. This property of CAN reflects the number of different graphs created during the simulation. This establishes all the graphs populated at least once.

2. Appearance. For any of the Ω X -items, there is a (0,1)matrix, i.e. X ω | X [i, j] where ω ∈ [1, ρ X ], that accounts for the appearance of X ω in X [i, j], such that:

(X ω | X) ij =    1 → X ij = X ω 0 → X ij = X ω
From this property of CMAS it is possible to identify when a given cluster is created.

3. Residence time distribution. The residence time of a given X ω is the amount of time that it appears continuously in a given trajectory; therefore, the residence time distribution (RTD) is the frequency distribution of the residence time for the whole ensemble of trajectories. It should be noted that the resolution of RTD is related to the value of τ employed. From this property of CAN it is possible to determine a lifetime of the desired graph if and only if τ is short enough.

Maximum-residence time distribution.

The maximum-residence time distribution (MRTD) is the frequency distribution of the maximum residence time value per trajectory. The schematic representation of the algorithm used to obtain RTD and MRTD is reported in Supporting Information, Table S1. Note that, only those values for which RT > 0 and MRT > 0 are useful for the analysis. From this property of CAN it is possible to determine the maximum time that a graph can be formed continuously.

5. Time-dependent probability and standard deviation. All the items in Ω X have associated 1D-matrices which account for the probability P ω [j] with the associated standard deviation σ ω [j] of the appearance of X ω in X [i, j] along the discrete time (j -1) τ for an ensemble of N-trajectories, such that:

P ω j = 1 N N ∑ i=1 (X ω | X) ij (3)
It corresponds to a reduction of dimensionality of X ω | X [i, j] by calculating the average over the trajectories. From this property of CMAS the population of a given cluster is obtained as a function of time .

6. Mean value or hierarchy. The average value of P ω [j] is expressed directly by:

P ω = 1 N • M N ∑ i=1 M ∑ j=1 (X ω | X) ij (4)
This value can be used as a criterion to identify whether an item is more important than others during the simulations. The item X µ is more relevant than X ν iff the P µ is greater than P ν . From this property of CAN it is possible to establish a hierarchy relative to a global population of a graph in the simulation.

7. Shannon's entropy. The disorder of X [i, j] can be measured by using Shannon's entropy formula, such that:

S X = - ρ X ∑ ω=1 P ω • ln P ω (5)
From this property of CAN it is possible to obtain a value that reflects the disorder of the system in the simulation.

8. Closure. As a matter of consequence from the property 6, the closure condition can be written as,

ρ X ∑ ω=1 P ω = 1 N • M ρ X ∑ ω=1 N ∑ i=1 M ∑ j=1 (X ω | X) ij = 1 (6)

IV. THREE-STATE KINETIC MODEL

We now consider the evolution of an ensemble of trajectories whose initial conditions were set in order to reproduce a specific physical condition. In the present case is the vibrational activation under microcanonical conditions. The following refers to the application of the general graph-theory approach described previously to the reactivity of an ensemble of microcanonical trajectories (all trajectories have the same internal energy). The method is general and can be applied without loosing generality to canonical simulations, for example. Since the system has at t = 0 enough energy to react, we now describe the different steps necessary to follow the graph modification in time and then obtain the kinetic properties.

The sample space Ω CAN covers all the physically possible graphs ρ CAN that can be formed considering the λ -atoms.

This number depends on the number of bonds of the λ -atoms and increases exponentially with their number. Nevertheless, during the dynamics only a reduced number of graphs are reached as a function of the total energy of the system, the starting activated graph (SAG), the total simulation time (Λ), and the number of propagated trajectories (N).

The selection of Ω CAN as a work-space avoids counting as different graphs those snapshots in which the system exchange identical functional groups (H, H + , CH 3 , etc). The features of Ω CAN can be extracted directly from the properties of the CAN [i, j] described previously iff the number of propagated trajectories is large enough. For example, the number of visited graphs Ω CAN (E ν , N, Λ) during the simulations for a given E ν , N and Λ is ρ CAN .

We should note that if we use the simple adjacency matrix space, Ω ADJ , the description of the dynamics could be problematic from a chemical point of view. In fact, two graphs representing the same chemical structure could have different adjacency matrices, and thus, some additional information should be included. The use of canonical label avoids this further complication.

A. Partitioning of the Work-Space and Definition of the States

The first element of an Ω CAN is the canonical label (CL) of the starting activated graph (SAG). It corresponds to the CL CAN 1 of the equilibrium starting geometry Q 0 . The visited graphs during the simulations can be divided into two regions. The first one covers all the connected graphs. Chemically speaking, it accounts for all the possible isomers or intermediate states before the system breaks eventually into two or more parts. The second zone encompasses all the disconnected graphs. In this region, the primary fragmentations happen. All the related information to this classification is stored in the CLU [i, j] array. For a given trajectory i at time (j -1) τ, the CLU =1 | CLU ij element is one when we have a connected graph and analogously CLU >1 | CLU ij element is one when a disconnected occurs, otherwise they are zero. The label "= 1" is used when one connected subgraph is detected while "> 1" when they are two or more. Figure 4 shows the resulting Venn diagram of Ω CAN accordingly to the partitioning we employ.

The primary fragmentation process of the SAG can be analyzed from the time-evolution of three different global states:

• Primary fragmentation state |PF . A trajectory is said to provide a primary fragmentation event if the system stays a given amount of time, ∆, in the disconnected region continuously. The ∆ value must be large enough to account for only those trajectories that once they reach the disconnected region never recross to the connected one. Of course, the time discretization, ∆/τ = β ∈ N, due to the storing frequency time holds (and in any case it will be valid since computer trajectories are always time discretized). The |PF assignment starts by building the primary fragmentation guess matrix (PFG), such as:

PFG [i, j] = CLU >1 | CLU [i, j] (7)
Subsequently, the (0,1)-matrix that represents |PF is obtained. A schematic representation of the adopted algorithm is reported in the Supporting Information, Table S2. Note that here we are dealing with gas phase fragmentation for which once the |PF is obtained it typically does not come back to the connected region. The approach can be in principle extended with physical conditions in which a chemical equilibrium is obtained between reactants and products with some (minor) modification. However, this will not be discussed in the present work.

Note that a trajectory to be classified as being in the primary fragmentation region has to be a disconnected graph for a given amount of time. So the PF state is not a set of disconnected graphs, but a condition that they must meet along with the time.

The use of ∆ implies a reduction in the dimensionality of the PF-matrix, N × (Mβ ), compared to PFGmatrix, N × M. This means that the trajectories can be analyzed only during Λ -∆ time. Then, a key point is the proper estimation of ∆. Using a too small value will increase the number of incorrect assignment of PFs, while a too large value will decrease the statistics in the analysis at the same computational cost. In other words, a larger value of ∆ will provide fewer ion-molecule complexes wrongly assigned as primary fragmentation.

To evaluate in the present case a reasonable value of ∆ and provide an approach which can be used for other problems, we consider the maximum residence time distribution of PFG [i, j], previously defined. This gives the information about the maximum time the system stays in the disconnected region before re-entering into the connected one. The MRTD of the PFG matrix for two simulations of CYC 00 and LIN 00 at 177 and 151 kcal/mol, respectively, is reported in the Supporting Information, Figure S1. Using ∆ = 5 ps the error is less than 1%, which can be considered as a good compromise and used in the following. Note that this allows to remove only 5 ps from the analysis, thus providing 15 and 35 ps that can be analyzed from trajectories long 20 and 40 ps, respectively. Using a larger ∆ will imply to run longer simulations or to have a smaller statistical sampling for the analysis (so a larger uncertainty).

• Three types of subsets can be distinguished. The first one is simply SP SAG = {CAN 1 } ⊆ SP, corresponding to the CL of SAG, with the corresponding appearance matrix:

SP SAG [i, j] = CAN 1 | CAN [i, j] (8)
This is the major contributor to the starting points state, nevertheless, during the microcanonical normal mode sampling of initial conditions, Q 0 + ∆Q, some other canonical labels are inherently populated and recognized as different in the analysis. Of course, they depend on the α parameter used in the construction of the connectivities. However, it will be simpler and more practical to keep α constant and elaborate a strategy to identify the other CANs. They can be of two kinds which we list in the following.

-SP C = {CAN 2 , CAN 3 , ..., CAN δ } ⊆ SP: ensemble of CL of connected (C) graphs that have a similar structure to the SAG. It includes those graphs barely populated directly from the SAG that contain unphysical atomic valencies and, in some cases, those graphs directly populated from the SAG by a fast dynamical process such as proton transfer -this is likely to occur when treating protonated peptides, for example, in which compact structures are common in the gas phase. The corresponding appearance matrix is:

SP C [i, j] = δ ∑ q=2 CAN q | CAN [i, j] (9)
-SP D = {CAN δ +1 , CAN δ +2 , ..., CAN η } ⊆ SP: set of canonical labels of disconnected (D) graphs formed as a result of bond breaking directly from SAG when E v increases. Typically, these disconnected graphs result from the X -H bond rupture due to the fixed value of α. A crucial requirement is that any of these graphs should not be involved in the primary fragmentation event, e.g., dehydrogenation or deprotonation. Therefore, the appearance matrix becomes:

SP D [i, k] = η ∑ q=δ +1 CAN q | CAN [i, k] • {J [i, k] -PF [i, k]} (10) where k ∈ [1, M -β ] and J [i, k]
is a all-ones matrix. The circle represents the Hadamard product between matrices. As can be seen directly, the SP D matrix shows a reduction in the dimensionality due to its dependence on PF.

The (0,1)-matrix of the starting point state can be directly obtained by summing all of its components and considering the reduction of the dimensionality:

SP [i, k] = SP SAG [i, k] + SP C [i, j] + SP D [i, k] (11)
The set of graphs chosen as the starting point state for CYC 00 and LIN 00 system are reported in the Supporting Information, Figures S2 andS3. As can be seen, a proper selection of SP C and SP D makes |SP mostly independent on the α value. Only the trajectories that after the vibrational normal-mode sampling used to generate initial conditions, provide (so for t = 0) a canonical label equals to SP were propagated, being SP [i, 1] = 1 ∀i. This procedure allows to control the activation process, avoiding that the automatic sis identifies some unphysical isomerization or fragmentation at t = 0. Finally, we can obtain the timeevolution probability and the associated standard deviation of |SP . An example of the time evolution of the different starting points defined previously (SP SAG , SP C and SP D ) for two simulations, CYC 00 at 217 kcal/mol and LIN 00 at 191 kcal/mol, is reported in Supporting Information, Figure S4. As expected, the SP SAG represents the major contribution to the starting-point state while SP C and SP D represent a limited fraction.

• Intermediate state |INT : all snapshots that are not classified neither |SP nor |PF . The (0,1)-matrix is defined as:

INT [i, k] = J [i, k] -SP [i, k] -PF [i, k] (12)
where k ∈ [1, Mβ ] and J [i, k] is an all-ones matrix.

As a matter of fact, the INT matrix also shows a reduction in the dimensionality due to its dependence on the PF-matrix.

There are two types of intermediate states, and notably:

i. Connected intermediates INT C : they are the set of connected graphs (isomers) reached by the system before it reacts, for which

INT C [i, j] = CLU =1 | CLU [i, j] -SP SAG [i, j] -SP C [i, j] (13)
or

INT C [i, j] = ξ ∑ q=η+1 CAN q | CAN [i, j] (14)
where {CAN η+1 , CAN η+2 , ..., CAN ξ } is the canonical labelling of the connected graphs reachable in the simulation. ii. Disconnected intermediates INT D : it is a disconnected area where the primary fragmentation requirement is not fully accomplished. These structures in fact will come back to a connected graph after a given time. In the case of ion fragmentation, this region is where ion-molecule complexes take place. More in general, this region identifies the formation of (transient) complexes, which can occur also in the case of neutral fragmentation and can be related to a roaming mechanism. 58,59 We thus have:

INT D [i, k] = CLU >1 | CLU [i, k] -PF [i, k] -SP D [i, k] (15)
We should note that in the present approach the same graph can be classified as INT or PF depending on its fate. For example, if inside the ∆ time a disconnected graph comes back to the connected region it is classified as INT, otherwise as PF. Given the definition of these three states it is possible to follow their evolution as a function of time from which rate constants can be extracted, as we will show in Section IV C. Before, we will show in the next section that also energetic information can be obtained.

B. Energy Decomposition

The energy information stored in the POT [i, j]-array can be decomposed into three different arrays according to the defined states, such as:

POT SP [i, j] = SP [i, j] • POT [i, j] POT INT [i, k] = INT [i, k] • POT [i, k] POT PF [i, k] = PF [i, k] • POT [i, k] (16)
The • denotes the element-wise product (Hadamard multiplication) as introduced before. Each term corresponds to the potential energy of the three states: starting point, intermediate and primary fragmentation, respectively. The same decomposition can be applied to KIN [i, j]-array. From these newly defined continuous arrays the energy gap between the states can be obtained.

The V -value can be used to provide an estimation of the relative stability of the two structures during the dynamics. The potential energy difference between CYC 00 and LIN 00 at 0 K is about 26 kcal/mol at RM1-D level of theory which is used in the trajectories, while from the simulations it is between 11 and 12 kcal/mole depending on the activation energy, as shown in Table I. This will correspond to an effective energy difference which will be discussed later when comparing the threshold energies obtained from the fragmentation simulations.

In Figure 5 we show for each state the average (over the different trajectories per each snapshot) potential energy of the three states as a function of time as well as the corresponding probability density function (PDF). As an example, we display values for simulations of CYC 00 and LIN 00 at two activation energies (177 and 151 kcal/mol, respectively). Note that by considering the energy difference between the minimum energy geometries of the two structures (26 kcal/mol), the two sets have the same excess energy..

In the CYC 00 system the three states are very-well resolved in energy. The maximum values of the PDF are located at 90.5, 103.5 and 121.5 kcal/mol for |SP , |INT and |PF , respectively. The LIN 00 has the |SP located 11 kcal/mole above the |SP of CYC 00 , while the |PF has the same energy as in CYC 00 . This reflects that the fragmentation products are in average very similar for the two systems. This is because LIN 00 is the most populated isomer during the fragmentation of CYC 00 (as we will detail later). The average potential energy of |INT in LIN 00 varies along the time because the system will isomerize going to the most stable graph, CYC 00 . That explains why the distribution goes below the |SP energy distribution. Similar behaviors were obtained for other internal energies.

The same energy decomposition can be applied to every canonical label. Thus, the mean values of the potential energy of a given graph CAN q can be determined by:

POT q [i, j] = CAN q | CAN [i, j] • POT [i, k] (17)
The same holds for kinetic energy. This approach will provide mean values of kinetic and potential energies of the states defined by the state partitioning. They can be used to determine the stability of the different states and/or canonical labels from an average and dynamical perspective.

C. Kinetic model for a three-state system By using the three states defined and determined as detailed previously, we can now set up a kinetic scheme from which "effective" microcanonical rate constants can be obtained. We call them effective because they reflect the simplification of a complex kinetic network to a simple three-state model. In other words, the intermediate and product states do not correspond to specific molecules or geometries, but to a global definition. The kinetic scheme can be constructed directly from the Venn diagram previously introduced. From it a three-state Markov model (shown in Figure 6) describes the statistical kinetic process of primary fragmentation in the gas phase.

The |SP state is represented by a filled blue circle to indicate that at t = 0, it is the only populated state. Previously we have highlighted that the choice of the ∆-parameter guarantees that once the system reacts, either from the initial state or from the intermediate state, it never returns back to reactant or intermediate state. This is a typical situation of gas-phase reaction, in which the fragmentations are not reversible since the probability that the two fragments encounter each other in vacuum is almost zero. This means that a given trajectory will cross the Γ 3 or Γ 2 dividing surfaces only once and there are no reverse rate constants from PF to SP or INT. As we have discussed in the definition of INT state, if two (or more) molecules are formed but they form back one molecule the state is labeled as INT. In other words PF state cannot by definition comes back to reactants, and the choice of ∆ is crucial. Of course, further fragmentations are possible, but they are not discussed here. They can be treated with the same approach using a specific fragmentation product as initial state.

The kinetic equations of this general scheme can be writ-ten as:

dN SP dt = k SI • N INT -(k IS + k PS ) • N SP dN INT dt = k IS • N SP -(k SI + k PI ) • N INT dN PF dt = k PS • N SP + k PI • N INT ( 18
)
where N SP , N INT and N PF denotes the abundance of the starting point, intermediate and primary fragmentation states, respectively while k SI , k IS , k PS and k PI are the corresponding "effective" rate constants which for the present simulations are microcanonical. In general, for each k ij , "i" indicates the destination state and "j" the source state. The diagonal elements can be expressed as:

k SS = -(k IS + k PS ) and k II = -(k SI + k PI ).
From now on, the word "effective" is suppressed from the discussion for simplicity. The solution of the system of linear first-order differential equations assuming that at t = 0 the only populated state is N SP (0) = N 0 , can be obtained using the Laplace-Carson transform method. [START_REF] Song | Theoretical Mass Spectrometry[END_REF][START_REF] Rodiguin | Consecutive Chemical Reactions[END_REF] The analytical solution is reported in the Supporting Information.

Since the energy is conserved and enough to pass the barrier, the limit for t

→ ∞ is N SP /N 0 = 0, N INT /N 0 = 0 and N PF /N 0 = 1.
The microcanonical rate constants are obtained by making a simultaneous fit using the solutions of the differential equations system from the population of the states obtained in simulations. The simultaneous fit was made using a home-made program developed in C++ through the modular scientific software ROOT. The MINUIT library was employed during the minimization procedure. The ROOT scripts are provided in a public repository (https://github.com/afperezmellor?tab=repositories).

To obtain all the kinetic parameters we employed two approaches: (A) all the rate constants are obtained from the simultaneous fit; (B) the primary fragmentation rate constant is obtained from the reaction flux, and the others from the fit. We now show how we obtained the rate constants using the two approaches and how they compare each other.

Full simultaneous fit (method A)

The microcanonical rate constants are obtained from each ensemble of trajectories as a function of the activation energy by fitting the populations of SP, INT and PF using the kinetic equations described above by a simultaneous fit. In Figure 7 we report the resulting fits together with the population evolution obtained from CYC 00 simulations at the different activation energies (the same from LIN 00 are reported in Supporting Information, Figure S5). As can be clearly seen, the fit is very good and thus the rate constants can be extracted safely. This shows that the simple three-state kinetic model can globally catch the fragmentation dynamics of this system. Figure 7. Populations of the three states (gray lines) at different vibrational energies (reported in kcal/mol) of CYC 00 system fitted by the simultaneous fit, method A. Blue: starting point state; black: intermediate state; red: primary fragmentation state.

Reaction flux and partial fit (method B)

In the second approach used (method B) we first estimate the most relevant reactive rate constant from the simulations by a direct measure using the flux. [START_REF] Grote | [END_REF] In fact, the primary fragmentation can be described as a reaction flux in sample space across the dividing surfaces, Γ 2 (∆) and Γ 3 (∆). The role of the ∆ parameter is crucial to minimize the flux through these surfaces. We should notice that Γ 2 (∆) and Γ 3 (∆) do not correspond to any saddle point geometry, but they represent abstract dividing surfaces. This has the advantage of getting rate constants which represent the overall kinetics, independently on geometrical details which can, in general, make impossible to conceive and analyze the full set of reaction pathways characteristic of fragmentation of complex systems. We will discuss how a specific pathway of particular interest can be studied out from the global pathway description in section V.

In this framework, the reactive rate constants can be ob-tained as a function of time from the flux of reactive graphs passing across the critical region divided by the remaining graphs in the departure state, such that:

k PI (t) = 1 τ |INT → PF |INT → INT k PS (t) = 1 τ |SP → PF |SP → SP (19)
In the systems studied here, there are only a few trajectories that react directly from the starting point, so the sampling will not be enough to get a converged value. Therefore, we have calculated only k PI directly from the flux.

The rate constants obtained via the flux method for CYC 00 simulations are reported in the Supporting Information (Figures S6 and S7). The rate constants are here time dependent and we observe a not negligible spread. At low energies this spread is larger, probably because the fragmentation is not fully converged. LIN 00 simulations show a similar behavior.

We have then used the position of the maximum value of the probability distribution functions as k PI and fitted the other rate constants. As previously, the fit results describe very well the simulation data (see Figure S8 in the Supporting Information). We should notice that this approach has the disadvantage that a huge number of trajectories is necessary to have a convergence in the flux calculations, while less are needed when using only the solution of the kinetic scheme from population evolution.

D. Rate Constants

The rate constants obtained at different energies and with the two methods are reported in Table S3. As we can notice, the two methods are in very good agreement providing very similar results.

For both systems, the rate constant values increase monotonically with internal energy. Only the k SI of CYC 00 are almost energy independent. In the CYC 00 system, the k IS is always larger than k PI . A different trend is observed for LIN 00 system, where at low activation energy k IS is higher than k PI and vice-versa for larger internal energy.

The k PI values for LIN 00 system are always higher than its counterpart CYC 00 at the same total energy. This observation is easily explained by the fact that the primary fragmentation state is found on average at the same potential energy level for both systems for the same total energy (see Figure 5). In addition, the starting point state CYC 00 is lower in energy by 11-12 kcal/mol (as obtained by the simulations at different activation energies) with respect to LIN 00 . Therefore, the energy gap between the two states, |PF -|SP , is smaller in LIN 00 simulations than in CYC 00 . We should remind that they share the same reactive basin and, then, the system LIN 00 will have more primary fragmentation events than its counterpart CYC 00 under the same total energy condition.

E. Threshold Energies

Using the two approaches described previously, we obtained the rate constants as a function of the excess energy for the two systems. The rate constants are obtained from microcanonical simulations and thus their energy dependence should follow the classical RRKM theory (also called RRK theory): [START_REF] Baer | Unimolecular reaction dynamics: theory and experiments[END_REF] k

(E v ) = ν 1 - E a E v s-1 (20)
where ν and E a are the effective frequency and threshold energy, respectively, while s is the number of vibrational degrees of freedom of the reactant. ν and E a are two adjustable parameters that can be obtained by fitting the rate constants at different energies via Equation 20. Figure 8 shows how rate constants behave as a function of energy, E v , with the corresponding RRK fits. As previously discussed, a clear energy dependence was observed for k IS and k PT in both CYC 00 and LIN 00 simulations, and k SI in LIN 00 simulations. Note that for k PS we have too few points to obtain a reasonable fit. As already noticed, the k SI in CYC 00 simulations do not show any energy dependence. This is probably due to the fact that the threshold energy is too low and that there is no time for energy randomization of the INT state before reaction.

From the RRK fits, we obtained frequency factors, ν, and threshold energies, E a , which are listed in Table IV.

As can be seen, the RRK model explains properly the evolution of the microcanonical rate constants, obtained from the simulations, as a function of the excess energy. As expected, the LIN 00 system needs less energy to form the primary fragmentation products than CYC 00 system. The difference between the threshold energies, 11 kcal/mol, corresponds to the energy difference between the starting point states obtained from the potential energy average (between 11 and 12 kcal/mol as a function of the activation energy). This result strengthens the picture for which the stability of the states is better described by the expectation value of the ) and energy thresholds (in kcal/mol) obtained following the methods to obtain the rate constants and then fitted using the RRK model.

potential energy from the simulations rather than from the difference between potential energy values at minima. Concluding, the proposed three-state model gives a general overview of the description of the primary fragmentation events providing quantitative information on rate constants and threshold energies.

F. Structures of the Intermediate States

The intermediate state is defined as a state which is different from the reactant but not (yet) primary fragmentation state. We now show that our approach can provide more molecular information on this state from graph-theory based analysis. First we can consider three subsets: linear connected intermediates, cyclic connected intermediates, and ion-molecule complexes (IMC). The first and the second subsets account for the linear and cyclic structures reached after the starting structure is activated. The latest are those disconnected structures that do not accomplished the primary fragmentation requirement.

The branching ratio of the intermediate state as a function of time as obtained at a given energy (the same total energy) from CYC 00 and LIN 00 simulations is shown in Figure 9. For the CYC 00 system, the linear structures represent almost 90% of the intermediate state while the cyclic structure is around 7% and IMC 3%. A completely different picture is observed for LIN 00 system. The linear structure drops down up to 47% while the cyclic formation grows up to 47%. The IMC represents around 7% of the intermediate state.

The explanation of such difference lies in the nature of the isomerization process. The CYC 00 system, which is the global minimum, has a cyclic structure. Once it is activated, the ring opens and thus the mostly formed intermediates are linear, such as LIN 00 . Note that in Figure 9 we show the ratio between three classes of intermediate states (INT) and thus the initial state (SP) is not included. After the ring opening, several proton transfers take place populating many linear structures. On the other hand, when the linear, LIN 00 , is activated, the main isomerization event is the formation of a cyclic structure, and notably CYC 00 , which is the most stable isomer.

It is possible to follow the time-abundance of the different structures of the intermediates. We report in Figure 9 the evolution of the most relevant ones (the corresponding struc- tures are reported in the Supporting Information, Figure S9, together with the other intermediates). We can notice that the linear structures mostly populated correspond to different tautomers: the proton transfer is a key reaction occurring to these intermediate structures, in particular for linear structures, following the mobile-proton model. [61][62][63] On the other hand, once a cyclic structure is formed this is much more stable than the linear counterpart. These processes are typical of peptide fragmentation and we have noticed some of them in previous simulations: 6,7,30,31,52,69,70,[79][80][81][82][83] here our approach is able to automatically characterize them and provide information on the different structures before their visual inspection.

The structures of the ion-molecule complexes (IMC) obtained during the simulations are reported in Figure S10 of the Supporting Information. They appear as pairwise fragments; one entity is charged while the other is neutral.

During the dynamic, they may interchange some functional groups, as observed here, resulting in different IMC structures: IMC 01 01 , IMC 01 02 , IMC 01 03 and IMC 01 04 . This is an indication that a roaming mechanism took place.

V. PROPERTIES OF A SPECIFIC PATHWAY

Up to now, the kinetic properties were obtained for the primary fragmentation state as a whole without any differentiation on their structure. We now show how the kinetic scheme can be adapted to obtain the desired properties for one specific path of interest. The approach is totally general, and of course if more specific products are interesting the procedure should be repeated for each pathway. Of course, one needs to have enough trajectories leading to a specific pathway to have statistically significative information.

When a primary fragmentation is identified, we can distinguish three classes of events:

1. Direct primary fragmentation : M + → M + 1 + M 2 . It occurs when the molecule reacts directly forming two entities of masses M 1 and M 2 . The information related to this events can be extracted from the MAS array. It is then useful to focus on a particular event, say M α = {M 1 , M 2 } which can be monitored during simulation using the properties described previously. The condition that a specific fragmentation occurs formally corresponds to :

MASD [i, k] = MAS α |MAS [i, k] • PF [i, k] (21)
Then, MASD (0,1)-array displays whether a given snapshot contains an event of direct primary fragmentation or not, and the time evolution can also be obtained.

Ion-Molecule formation : M

+ → M + 3 • • • M 4 + M 2 where M + 3 + M 4 = M + 1 .
This process occurs when a complex is obtained during the fragmentation that gives M + 1 . This event can be followed through the CMAS array which, as previously seen, takes into account the formation of clusters. Similar to what done for direct events, and considering the same M α = {M 1 , M 2 } event, it can be stated that:

MASC [i, k] = MAS α |CMAS [i, k] • PF [i, k] (22)
Note that this is not an intermediate ion-molecule complex, but a case in which one (at least) product is itself a complex.

3. Higher-order of fragmentation event: they are structures which are neither direct fragmentation nor ionmolecule formation events. These events are predominant when the activation energy is too high and the system breaks suddenly in many pieces.

The formation of a given mass channel M α = {M 1 , M 2 } via direct primary fragmentation can be obtained by:

MASG [i, k] = MASD [i, k] + MASC [i, k] (23)
Figure 10. Markov diagram for a particular four-state system that explains the formation of a given fragmentation product by primary fragmentation in the gas phase.

It should be noted that M 2 can fragment subsequently into different pieces depending on its nature and excess energy.

In that case, one can follow simply M + 1 and M + 3 • • • M 4 as primary fragmentation products, disregarding the fate of the neutrals. At this end one should use the information on the fragmentation products. These events were not found here due to the low mass of the neutrals and relatively low fragmentation energy. This approach will provide the different fragmentation channels, in terms of clusters and fragments (an example is given in Table S4 of the Supporting Information). This information will be useful to build the theoretical mass spectrum, as we will show in section VI. Before, we will show how the kinetic parameters of a specific fragmentation product can be obtained.

Once we monitor the appearance of specific fragmentation channels, it is then possible to determine their unimolecular kinetics. Similar to how it was done to explain the global process of primary fragmentation, to extract the rate constants we first built a kinetic model with four states, shown in Figure 10. The primary fragmentation state is now divided into two new states, |A and |B . The first takes into account the process of interest, for example, the formation M α = {M 1 , M 2 } while the second all the other processes, so that |B = |PF -|A .

After establishing all possible interconnections between the states and considering that once the system reacts it does not return to the precursor state, the system of ordinary dif-ferential equations will be :

dN SP dt = k SI • N INT -(k IS + k AS + k BS ) • N SP dN INT dt = k IS • N SP -(k SI + k AI + k BI ) • N INT dN A dt = k AS • N SP + k AI • N INT -k BA • N A dN B dt = k BS • N SP + k BA • N A + k BI • N INT (24)
The analytical solution (reported in the Supporting Information) of this system was obtained in a similar way to that of three-state model using the Laplace-Carson transform method. It is important to note that there are four totally independent populations N SP , N INT , N A , and N B and seven rate constants to be determinated: k IS , k SI , k AS , k BS , k AI , k BI and k BA . However, two have already been obtained in the three-state system k IS and k SI together with the sum of two of them as, k PS = k AS + k BS and k PI = k AI + k BI . Therefore the problem is reduced to four populations and three rate constants, k AS , k AI and k BA . The rate constant k AI accounts for the formation of the mass channel through a global intermediate while k AS does the same but from the starting point. Finally, the rate constant k BA corresponds to the secondary fragmentation of the channel of interest if there is any.

A simultaneous fit of the populations, following the method A, was performed to obtain the rate constants involving the |A state, k AS , k AI and k BA . The remaining rates constants, k IS , k SI , k PI and k PS , are constrained to the values obtained during the simultaneous fit, by using the approaches A or B, of the three-state scheme.

We applied this approach to two primary fragmentation products, which were the most abundant observed: (i) the loss of neutral CO and (ii) the formation of methaniminium cation (H 2 NCH + 2 ). As previously, we obtained the rate constants from the time evolution of the states and then we fit their energy dependence using the RRK model. The RRK fit results are reported in Table S5. As for the global fit, the values obtained with method A and B are very similar. For simplicity we will comment in the following the values obtained with method A.

Interestingly, in the case of CO-loss we can follow two reactions: (i) the formation of the CO-loss product from an intermediate state (characterized by ν AI and E a AI ) and (ii) the degradation of it forming other products (characterized by ν BA and E a BA ). For the |INT → |A reaction, the threshold energy (E a AI ) is slightly lower than the corresponding global threshold energy (E a PI ), which is what expected since this channel is the most abundant one. The same picture is obtained for CYC 00 and LIN 00 simulations, just the barrier of these lasts is lower. Notably, the differences in threshold energy reflect also in this case the difference of 11-12 kcal/mol obtained from average potential energy analysis.

For the H 2 NCH + 2 formation we do not observe any secondary fragmentation and the energy threshold from CYC 00 simulations are about 14.5 kcal/mol higher than the one obtained from LIN 00 simulations. This is slightly higher than the average energy difference between the two states, but not large enough to suggest that other processes are relevant.

VI. THE THEORETICAL MASS SPECTRUM

We now show how this approach can finally provide a theoretical mass spectrum. More detailed description of the reactivity will be left to specific studies as done previously for other systems. 65,66 Experimentally this is obtained by counting the number of structures as a function of their mass-tocharge (m/z) ratio. Theoretically we can obtain it in the same way using the rigorous definition of a primary fragmentation product given before. Also an estimation can be made by counting the elements in the Ω CMAS work-space. Every item in Ω CMAS includes the set of cluster mass formed during the simulations.

All the different items ρ CMAS are directly related to the possible reactions. For any of these different items ω ∈ [1, ρ CMAS ], the charge distribution is calculated a posteriori from the expectation values of the electrostatic potential (ESP) of the molecule on a uniform distribution of points as implemented in Mopac. [START_REF] Stewart | MOPAC, version 5.022mn based on MOPAC 5[END_REF] Experimentally, the most abundant m/z products observed in fragmentation of protonated cyclo-di-glycine are: 60 90.0, 87.0, 70.1, 59.0, 58.1, 41.9 and 30.1. We obtained all of them from the different simulations (results for the two initial species and the different excess energies are summarized in Tables S6 andS7 of the Supporting Information) as well as other products. The description of the observed ions is reported in Table S8. Note that from simulations we obtain a high-resolution spectrum with no additional cost. Furthermore, it would be eventually possible to consider isotopic abundances, assuming that the fragmentation process is not largely affected by difference in atomic masses (which is what generally is done when doing isotopic labeling in mass spectrometry).

Similarly to experiments, it is possible to report the relative abundance of the species of interest as a function of activation energy (an example is reported in the Supporting Information, Figure S11).

VII. CONCLUSIONS AND OUTLOOKS

In this work, we have reported a new application of graph theory to analyze an ensemble of reactive trajectories of a complex system. We focus our attention on the gas-phase unimolecular dissociation related to mass spectrometry, but the approach is totally general. One key aspect is the use of the canonical labeling which allows to automatically classify as the same molecular structures two graphs with different adjacency matrices. A flowchart summarizing the key steps to obtain the main information to characterize the reactivity using graph-theory and kinetic analysis is reported in Figure 11.

We should note that recently the spectral graph-theory was used to analyze molecular reactivity, 50,84 which could be less time consuming for much larger systems. Presently the We report the different matrices obtained in their sequence and the different parameters, states and properties as defined in the text.

graph-theory analysis using the canonical labeling is much faster than the production of trajectories, and thus the computational time for the trajectory analysis was not a problem. However, when using spectral-graph theory, as in the SPRINT coordinate approach 45 which considers the permutation of identical atoms, a matrix diagonalization is needed which would also be time consuming for very large systems. Future comparisons of the different approaches will surely be useful given the increasing interest of applying graph theory to study chemical reactivity.

Based on the present approach, it is possible to obtain kinetic properties, in terms of rate constants and, after a typical RRK fit, threshold energies. Here we assumed that the products cannot form back the reactants, as it is pertinent in gas phase fragmentation, but a particular care is given in identify transient species which can be obtained, as the ionmolecule complexes, this being relevant to point out possible roaming mechanisms. The method, however, could be extended to physical conditions (and typically reactions in solution) in which a dynamical equilibrium between reactants and products is reached. The main limitations will be in having enough direct dynamics trajectories to have a good statistical sampling. Also the kinetic scheme should be adapted in this case. Future applications and developments of the present method in that direction are surely welcome.

Both general three-state and specific four-state models can be used, this last allowing a more specific analysis of a particular fragmentation product. We should also note that we obtained rate constants also from the flux: results are similar to the simple fit of the population while this last needs less trajectories to have a good convergence. We thus suggest in future to use it when possible.

A particularity of the approach is the general definition of an intermediate state, which is rather abstract (this is useful to obtain rate constants of fragmentation products) but molecular details can also be extracted. Specifically to ion fragmentation is the possibility to identify ion-molecule complexes which are often suggested when explaining products observed in mass spectrometry.

Finally, we have coupled the graph theory analysis to a charge localization method leading to a full theoretical mass spectrum which can be obtained automatically. We have applied it to a prototypical system representative of peptide fragmentation, we now plan to use it in more complex and new systems in order to describe not only qualitatively but also quantitatively (thanks to the access to threshold energies) the fragmentation processes.

Figure 1 .

 1 Figure 1. Schematic representation of the two structures used in chemical dynamics simulations. CYC is for the cyclo-di-glycine and LIN for its linear isomer. Both are protonated forms. The label 00 means that they are the most stable structures. The charge is arbitrarily placed on the C atoms, another valence structure can be drawn with the charge on the O atoms and C=O double bonds. The simulations are done in the molecular orbital framework so this will not have any impact.

Figure 2 .

 2 Figure 2. A snapshot taken from the simulations used to illustrate the graph decomposition method used in the text. CL stands for the canonical label, while dclus represent the cut-off distance used to identify the cluster formation; see text for details.

  Once all connected subgraphs are identified in ADJ ij , their corresponding vertices are sum-weighted by the atom mass. Then, each MAS ij -element corresponds to the set of mass recognized in ADJ ij . For example, MAS ij = {M1, M2, M3} where M1 > M2 > M3 if three different connected subgraphs were detected in ADJ ij , as in Figure2. In the present work we used the following masses: H C N O M (uma) 1.00782 12.00000 14.00307 15.99490 which correspond to what used in the dynamics.

  The cyclomatic number of the graphs can give the number of rings present in the structure, c = |E| -|V| + D where |E|, |V| and D are the number of edges, the number of vertices and the number of disconnected components

Figure 3 .

 3 Figure3. The PDF of the distance between disconnected graphs normalized at its maximum value as obtained from two simulations: CYC 00 with an internal energy of 177 kcal/mol (blue line) and LIN 00 at 151 kcal/mol (green line). The two vertical lines highlight the maximum of the PDF and the dclus value used. The sketch summarizes the origin of the major contributors to the peak of the PDF.

  Information Matrix The position of the atoms XYZ ij The potential energy of the graph POT ij The kinetic energy of the graph KIN ij The adjacency matrix of the graph ADJ ij The total energy of the graph TOE ij The partition of the graph π The canonical label of the graph CAN ij The number of connected subgraphs in the graph CLU ij The adjacency matrices of the subgraphs The canonical labels of the subgraphs The distances between connected subgraphs The masses of the connected subgraphs MAS ij The masses of the clusters CMAS ij The number of cycles contained in the graph CYC ij

Figure 4 .

 4 Figure 4. Venn diagram of Ω CAN displaying the definition of the states. Γ establishes the boundaries between different states.

  Starting point state, |SP . It corresponds to the set of CLs which have the same chemical structure of the geometry considered as initial state: SP = {CAN 1 , CAN 2 , ..., CAN η }. This problem can arise because the initial structure is activated randomly and differently for each trajectory. The automatic algorithm can find a different connectivity and thus a different CAN, while they represent the same physical state and one should analyze them together.

Figure 5 .

 5 Figure 5. Average potential energy as the function of time (left panel) and corresponding probability density function (right panel) as obtained from chemical dynamics simulations at two activation energies : a) CYC 00 at 177 kcal/mol and b) LIN 00 at 151 kcal/mol.

Figure 6 .

 6 Figure 6. Markov diagram for the three-state system that models the unimolecular dissociation in the gas phase.

Figure 8 .

 8 Figure 8. Microcanonical rate constant values (symbols) of the different processes as a function of the internal energies. The colored line represents the RRK fit from which the energy threshold and the effective frequency factors are obtained.

Figure 9 .

 9 Figure 9. Branching ratios : Upper panels, decomposition of the intermediate state between linear, cyclic, and IMC structures at the same total energy. a) CYC 00 simulations at 177 kcal/mol. b) LIN 00 simulations at 151 kcal/mol; Lower panels, the three most relevant structures of the intermediate state as obtained from simulations: a) CYC 00 at 177 kcal/mol and b) LIN 00 at 151 kcal/mol.

Figure 11 .

 11 Figure 11. Schematic flowchart summarizing the present approach. We report the different matrices obtained in their sequence and the different parameters, states and properties as defined in the text.

Table IV .

 IV IS 0.110 ± 0.014 0.111 ± 0.016 0.005 ± 0.0006 0.005 ± 0.0009 E a IS 30.480 ± 0.517 30.518 ± 0.561 13.845 ± 0.561 13.725 ± 0.675 Effective frequencies factors (in fs -1

	CYC 00		LIN 00
	Method		Method
	A	B	A	B
	ν ν SI E a SI		0.001 ± 0.0004 0.001 ± 0.0007 11.147 ± 1.116 11.850 ± 1.619
	0.263 ± 0.037 PI 37.859 ± 0.919 36.169 ± 0.919 26.375 ± 0.919 26.915 ± 0.362 0.113 ± 0.009 0.133 ± 0.013 ν PI 0.400 ± 0.096 E a

ACKNOWLEDGMENTS

This work was supported by ANR DynBioReact (Grant No. ANR-14-CE06-0029-01). We thank Dr. Anne Zehnacker, Dr. Alejandro Gutierrez Quintanilla, Dr. Latevi Max Lawson Daku, Dr. Arnulf Rosspeintner, and Prof. Thomas Burgi for useful discussions and suggestions on the manuscript. We also would like to thank the late Prof. William L. (Bill) Hase for several discussions we will clearly miss in the future.

SUPPLEMENTARY MATERIAL

Supplementary material is available, where we report: (I) eleven additional figures, (II) eight additional tables and (III-IV) the analytical solutions of 3-state and 4-state kinetic systems.

The scripts used in the kinetic fits are reported in a public GitHub repository (https://github.com/afperezmellor?tab=repositories).

AUTHOR'S CONTRIBUTIONS

All authors contributed equally to this work.

DATA AVAILABILITY

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