Theoretical Calculations of the Multistep Reaction Mechanism Involved in Asparagine Pyrolysis Supported by Degree of Rate Control and Thermodynamic Control Analyses

: A computational study on the mechanisms of reaction for the pyrolysis of asparagine is presented. A density functional theory (DFT) study at the ω B97XD / 6-311G(d,p) level was performed to analyze the di ﬀ erences in two reaction mechanisms: (i) the formation of ﬁve-membered cyclic products: maleimide and succinimide, and (ii) the more classical, six-membered cyclic products (diketopiperazine species) which are common in the pyrolysis of many other amino acids. The e ﬀ ect of temperature was included in the calculations at 300 ◦ C or 625 ◦ C, as required. Moreover, a detailed study based on the degree of rate control and thermodynamic control of the proposed mechanism for the formation of maleimide and succinimide is also presented. Results show that, for asparagine, the ﬁve-membered ring formation is the preferred process instead of the six-membered cycle (32 kJ / mol of Gibbs free energy di ﬀ erence between them at the ﬁrst cyclization step); therefore, the polymerization is favored. On the other hand, the rupture of the polymer represents the highest energetic barrier ( ∆ G ‡ = 281 kJ / mol) and the most inﬂuential process in the overall rate of the reaction. These results are in good agreement with the experimental evidence.


Introduction
Pyrolysis is defined as the chemical reactions occurring exclusively by the input of thermal energy to a reactant during heating processes above certain temperatures [1]. There is a wide range of applications for the pyrolysis techniques, which are normally grouped in the following two approaches: (i) applied pyrolysis, which is more concerned with the yields of pyrolysates [2] and has many large-scale industrial uses, such as the synthesis of raw materials from the petroleum refining industry [1], and the pyrolysis of biomass to produce fuels and chemicals [3]; and (ii) analytical pyrolysis, which is applied to the characterization of a particular analyte, the study of thermal reactions of simple organic molecules, or as a basis for mechanistic interpretations [2,[4][5][6][7]. For instance, pyrolysis-gas chromatography/mass

The Formation of Maleimide and Succinimide Products: Mechanism of Polymerization-Rupture
Maleimide and succinimide, five-membered-heterocyclic products, have been experimentally obtained as main products during pyrolysis of asparagine, after the formation of a polysuccinimide structure [8][9][10]12]. The mechanism for the process is proposed in this work, as shown in Scheme 1.
The proposed mechanism is separated in two main processes in concordance with the experimental observations: the polymerization and the rupture of the polymer. The model studied in this work regards the synthesis of an oligomer with four monomers of asparagine, obtaining an intermediate (INT12) with three succinimide-like rings in its structure; from this intermediate, three subsequent individual ruptures finally lead to the products (P1 and P2). Each addition of monomer completes after four elementary steps: a condensation between carbonyl and amine groups, subsequent dehydration of the glycol group formed, the intramolecular formation of a five-membered ring, and the elimination of a molecule of ammonia. On the other hand, each rupture generally achieves with two elementary steps: the removal of a residue and a subsequent tautomerization.

Optimization and Thermodynamic Properties
The structures of all transition states (TS) and all stationary points (including reactive R, intermediates INT, products P, byproducts H 2 O and NH 3 , and residues, Res) were optimized, and frequency calculations were performed to obtain the corresponding thermal properties of each structure. Finally, the program Goodvibes v2.0.3 [14] was used to obtain the quasi-harmonic corrected Gibbs free energies presented in Table 1 and plotted in Figure 1. An average difference of 14 kJ/mol was obtained when comparing to the non-corrected Gibbs free energies obtained by default in Gaussian16. Scheme 1. Condensed mechanism of reaction proposed for the formation of maleimide and succinimide during the pyrolysis of asparagine. Table 1. Thermodynamic parameters for the mechanism of polymerization-rupture of asparagine obtained through DFT calculations: wb97xd/6-311g(d,p). A quasi-harmonic correction was additionally applied using the program Goodvibes v2.0.3.  The colored bands (white and grey) separate the three additions of monomer during polymerization (with four elementary steps for each one), and the three individual ruptures during the entire rupture process of the polymer (with two elementary steps each one). Curved arrows represent formation of by-products, residues, or products that come out of the profile. Some labels are not shown in the figure due to clarity reasons.

Degree of Rate Control and Thermodynamic Control
In order to identify the rate-determinant process of this multistep profile ( Figure 1) for formation of the maleimide and succinimide, which is difficult to discern by simple observation of the pattern, the method proposed by Uhe et. al. [15] was used; thus, degrees of rate control (X RC ) for transition states and degrees of thermodynamic rate control (X TRC ) for minimum energy configurations (reactive and intermediates) can be calculated using Equations (2) and (3), respectively. The matrix of exponential terms (M jk ) obtained from the corresponding Gibbs free energies is shown, in part, in Equation (1) (See the complete M jk matrix in Table S2 in the supporting information), while the corresponding results for degrees of rate controls are shown in Table 2.
The degrees of rate constants indicate that TS13 is the most rate-influent transition state in the entire profile, with a weight of 99% in the global reaction rate (X RC,13 = 0.999). This transition state represents the starting point of the rupture of the oligomer, which is also the most demanding step from the mechanistic point of view (see Scheme 1), showing two bond breakings and a relatively large distance for proton transfer. Additionally, the immediately preceding intermediate of the mentioned rupture (INT12) accounts for almost 100% of the total thermodynamic rate control of the entire reaction (X TRC,13 = 0.971), and expresses the fact that the beginning of the rupture is definitely the most important rate-influencing process.

Reaction Force Analysis
Following the electronic energy over the reaction coordinate from INT12 to INT13 (step 13) by IRC calculations, an interesting behavior was obtained (see Figure 2), whereby the relaxation of the molecule after achieving TS13 showed a kind of plateau in the middle of the curve, which caused the existence of additional inflection points. The reaction force approach was applied, and six different works (w j ) were identified in the reaction force plot. Works w 1 and w 2 are located before TS13; they represent the 66% and 34% of the energy needed to achieve the transition state, and can be easily related to geometrical and electronic rearrangements, respectively. Works w 3 and w 4 , located after TS13, are also related to the electronic and geometrical reorganization in the product formation, respectively. Finally, works w 5 and w 6 located around the part corresponding to the last inflection point can be associated with a barrierless process where the proton transfer occurs from nitrogen to the oxygen. As mentioned earlier, the first rupture process represents a difficult step from the mechanistic point of view. Looking at the structures involved during the rupture process through the reaction coordinate (see Figure 3), one can see that the proposed two bond breakings occur before the transition state, breaking the ring and separating the residue from the molecular matrix, while the proton transfer apparently happens on the final stage, just in the part of the reaction force plot where there is a vertical asymptote (around IRC = 0.9).
To get more insight into this rupture, the geometrical changes through the intrinsic coordinate of this elementary step were studied (see Figure 4). The bonds that break (blue and green lines) start to lengthen at almost the same time, just before entering to the second region of the plot (associated to work w 2 ); then, at the beginning of the fourth region, the bonds are totally broken, although the distance between those atoms slightly increases. The hydrogen that is transferred keeps joined to the nitrogen at a constant range (red line) until almost the end of the fifth region of the plot; the hydrogen transfer then occurs, which is noticeable because the bond length of N-H rapidly grows at that point. Then, the distance between the oxygen and the hydrogen that will form a bond (yellow line) shows a very interesting plot: it decreases in the first region until a minimum, after which a slight increase occurs when the rupture of the ring starts; then, it continues falling until the proton transfer at the end of the fifth region occurs.  When the distance O-H (yellow line) decreases to approximately 2.5 angstroms, the tension in the ring is strong enough to break the bond N-C (green line), and the bond C-C (blue line) is simultaneously broken in order to redistribute the electronic populations on the separated molecular group. The increment in the distance O-H at that point assumed to be an effect of the molecular movements during the ring ruptures. The tensioned ring mentioned is stabilized by the hydrogen bond interactions between the oxygen and hydrogen atoms.

NBO Analysis
When the changes in the NBO atom charges-obtained from the respective natural bond orbital (NBO) calculations-were studied from INT12 to INT13 (see Figure 5), it was possible to identify that the larger charge transfers occurred around the transition state (IRC ≈ 0.45), while in the rest of the intrinsic coordinate, the NBO atomic charge kept almost constant in the selected atoms. Relatively large changes of charge also occurred around the vertical asymptote of the reaction force plot (IRC ≈ 0.9). In the first section of the plot (associated with work w1), there are no significant changes of charge; therefore, this work can be related to geometrical rearrangements, as mentioned earlier. The second and third sections show substantial charge changes, representing the rupture of the ring: N(2) and C(3) are the atoms that lost more electronic density, while O(9) and C(7) are the ones that gained more negative charge; these four atoms are located in the extremes of the double rupture (bonds C(7)-N(6) and C(4)-C(3) are broken), while atoms C(4) and N(6), located between the two broken bonds, are practically not affected by the rupture (in terms of atomic charge), which means that their charges are immediately compensated, and therefore, it can be assumes that both bond breakings occur practically at the same time.
In the fourth and five sections, charge changes are practically negligible except for the small area around the vertical asymptote, which represents the final hydrogen transfer: O(9), that takes the hydrogen and loses electronic density due to the bond formation; and N(2) that drops the hydrogen and gains electronic density due to the rupture of the bond. It is interesting to notice that C(7) and C(3) also relaxed, losing and gaining electronic density, respectively; thus, the evidence suggest that they were sharing the respective positive and negative charge densities before the proton transfer. Sharma et al. [12] generalized the formation of 2,5-diketopiperazine (DKP) products based on the pyrolysis reaction of some amino acids, such as proline and glycine, describing a reaction with two main processes: a dimerization and a cyclization [12]. The mechanism involved in this process was evaluated recently in detail by us [16]. The DKP mechanism (DKP mech) consists of a dimerization in the first stage, which is reached after the intermolecular condensation between a carbonyl group and an amino group, and is followed by a dehydration of the glycol group just formed. The cyclization process occurs by similar elementary steps, but the condensation step is intramolecular, and the dehydration from the glycol group created leads to the 2,5-diketopiperazine product (Scheme 2).

Scheme 2.
Mechanism of reaction proposed for formation of 2,5-diketopiperazine during pyrolysis of asparagine.
According to the steps just mentioned, a Gibbs free energy profile was built by the theoretical determination of thermal properties (DFT optimization and frequency calculations) at 300 • C for all transition states (TS) and stationary points (reactive, R; intermediates, INT; and products, P) involved in the mechanism; a quasi-harmonic correction was additionally applied using the Goodvibes v 2.0.3 program [14] (an average difference of 9 kJ/mol was obtained when comparing with non-corrected Gibbs free energies). The energy profile is shown in Figure 6. Since the change in the Gibbs free energy for the global reaction is negative (∆G reaction = −52 kJ/mol), it can be stated that the reaction is spontaneous. The dimerization process contains the highest energetic barrier in the entire profile (the barrier from 2R to TS1), and according to the degrees of rate control and thermodynamic rate control calculated (See Table 3), this simple step is stablished as the rate-determinant process of the entire reaction (See the corresponding M jk matrix in Table S1 in the supporting information). Therefore, the amino acid is required to overcome a barrier of 272 kJ/mol corresponding to the first step of the dimerization (the condensation step), while the subsequent elimination of water is considered a more straightforward step. The first difference between the two profiles described in Figures 1 and 6 is noticed in the cyclization process, especially from INT2 to INT3 in both cases. Looking at the respective mechanisms (Schemes 1 and 2), it is interesting to notice that the first and second steps are the same for both reactions, and the primary difference between them is the formation of a five-membered ring (∆G ‡ = 171 kJ/mol) in the mechanism of polymerization-rupture instead of the six-membered ring (∆G ‡ = 203 kJ/mol) obtained in the DKP mechanism, being the polymerization process kinetically favored by about 32 kJ/mol (Scheme 3). The difference in the Gibbs free energy activation barriers would constitute an excellent answer to the question about why, during pyrolysis of asparagine, the 2,5-diketopiperazine product is actually not obtained, but maleimide and succinimide molecules are obtained instead, as reported in experimental works [8][9][10]12]. Since the atoms involved in each process are not the same, there are evident differences from the mechanistic point of view (see Scheme 3). To obtain the 2,5-diketopiperazone product, an amine group reacts with a carboxylic group, while for the mechanism of polymerization-rupture, two amide groups are reacting with each other. Therefore, they are two different reactions. Scheme 3. Atoms of the dimer of asparagine that are involved in the formation of a six-membered ring (blue) according to mechanism for formation of 2,5-diketopiperazine (DKP), and a five-membered ring (red) according to the proposed mechanism of polymerization-rupture of the polymer (Polym).

Reaction Force Analysis on the Cyclization Process
The total electronic energy through the reaction coordinate of these two reactions was obtained by the respective IRC calculations (see Figure 7). The reaction force approach and determinations of the corresponding works were done in order to analyze the differences between these two cyclization processes. At first glance, the shape of the electronic energy plot corresponding to the mechanism of formation of 2,5-diketopiperazine is different from the plot corresponding to the mechanism by polymerization-rupture. This fact is in accordance with the already-mentioned differences from the mechanistic point of view. The electronic activation energy needed for asparagine to obtain the 2,5-diketopiperazine product (∆E ‡ = 175 kJ/mol) is practically the same as the barrier that corresponds to the polymerization-rupture mechanism (∆E ‡ = 174 kJ/mol), which is proposed for the production of maleimide and succinimide (the experimentally identified products of pyrolysis of asparagine). Then, the differences between the free energy commented on before are mainly associated with changes in the activation entropy of these two processes (∆ T∆S ‡ = 32 kJ/mol), which points out that the process involved in the polymerization-rupture mechanism implies a significantly greater loss in rotational and vibrational degrees of freedom than DKP formation. Figure 7. Total electronic energy, reaction force plots, and works, through the intrinsic reaction coordinate of the corresponding cyclization processes for asparagine: by the mechanism of polymerization-rupture (red), and by the mechanism of formation of 2,5-diketopiperazine (blue).
Despite the similarity between the energetic barriers, when values of work are analyzed, it is interesting to note that the mechanism by polymerization-rupture achieves the present cyclization step with very similar contributions of geometrical and electronic rearrangements (ratio w 1 /w 2 : 1.2), while the DKP mechanism has a higher contribution of geometrical rearrangements (ratio w 1 /w 2 : 8.7). That reasoning supports the explanation of why asparagine prefers to form the five-membered ring instead of the six-membered one.
Finally, the overall reaction free energies for both reactions are negative: −52 kJ/mol for the DKP formation and −87 kJ/mol for the polymerization-rupture reaction, which suggests that the final five-membered products (maleimide and succinimide) are more stable than the corresponding DKP product.

Optimization of Structures and Determination of Thermodynamic Properties
Density functional theory (DFT) has been extensively used by the scientific community to perform theoretical calculations on diverse molecular models and reactions [17]. Here, the Gaussian 16 suite of programs [18] was used to calculate the optimized geometry configurations and frequencies of all reactants, intermediates, transition states, and products involved in each of the mechanisms of the reaction we studied. Since reactions include non-covalent and nonbonded interactions between atoms, the long-range dispersion-corrected DFT functional ω B97XD [19] was adopted as the level of theory along with the 6-311G(d,p) basis set [20,21], whose polarized functions are expected to accurately describe the electronic configuration of atoms in the reacting zones. The validity of the obtained molecular arrangements was verified by analyzing the eigenvalues in the resultant Hessian matrix: minimum energy configurations were selected when no imaginary eigenvalues were obtained, while transition states had to contain a unique imaginary eigenvalue in the entire matrix (the transition vector corresponding to this imaginary frequency also had to be coherent to the expected reaction's elementary step).
A temperature of 573.15 K was set up in all frequency calculations to obtain thermodynamic parameters adequately-related to reported pyrolysis experiments reported in [12]. Absolute Gibbs free energies for each stationary point were obtained with the Goodvibes v 2.0.3 program (non-default settings: −t 573.15; −q grimme), a python program that recalculates all partition functions and applies a quasi-harmonic approximation to the vibrational entropy [14]; in this case, the free-rotor approximation, as proposed by Grimme [22].

Determination of Degrees of Rate Control and the Rate-Determinant Process
From the corresponding Gibbs free energy profiles, degrees of rate control for transition states, and degrees of thermodynamic rate control for minimum energy configurations were calculated according to the procedure described by Uhe et al. [15] (with some modifications in only notation). For a generalized sequence of n steps, Equations (2) and (3) were used to calculate the degree of rate control of the i-th transition state (X RC,i ) and the degree of thermodynamic rate control of the i-th minimum energy configuration (X TRC,i ), respectively: (2) where M represents the exponential term according to each combination of sub-indexes; G 0 TS is the Gibbs free energy of each corresponding transition state in the profile; G 0 Min is the Gibbs free energy of each corresponding minimum energy configuration (including the reactant and excluding the product); R is the gas constant; and T the temperature of the process. δG is either the difference of Gibbs free energy of the overall reaction (∆G r ), or zero, depending on the corresponding sub-indexes j or k, as specified in Equation (4):

Reaction Force Analysis
For the elementary steps of interest, the intrinsic reaction coordinate (IRC) algorithm [23] was used at the same level of theory. The total electronic energy (E) over the normalized reaction coordinate (ξ) was determined from the forward and reverse IRC profile, starting from the rate-controlling transition state optimized structure. After that, the reaction force (dE/dξ = −F) was calculated from the numerical derivative of the energy with respect to the reaction, which was used to build the reaction force profile corresponding to the rate-controlling process of the selected mechanism [16,[24][25][26][27].
In the reaction force profile, the reaction coordinates' points correspond to a cut across the horizontal axis (F = 0), and their maximum or minimum leads to the determination of regions (j-th) where the energy of the process can be analyzed in separated-way using the associated works (w j ). Values for works were calculated by numerical integration, as the area under each region (from critical point "a" to critical point "b": w j = − b a Fdξ). Those works were then used to determine how the energy of the process decompose into geometrical or electronic rearrangements, which relates to each region of the F profile.
Additionally, from IRC calculation outputs (that contain the information of all atomic coordinates in each point of the IRC profile), geometrical parameters in this case, some selected bond lengths between atoms of interest were followed through the reaction coordinate leading to additional analysis.

NBO Analysis
In order to gain more insights on the elementary steps of interest, NBO calculations [28] were performed for a sequence of different geometry configurations related to different intrinsic reaction coordinates (IRC points = A, B, C, ...). NBO charges (Q atom IRCpoint ) of selected atoms in the different IRC points were obtained, and comparing them against the NBO atomic charges corresponding to their immediately preceding IRC point, allowed the determination of consecutive changes of charge (δQ atom ∆IRC , where ∆IRC = B -A, C -B, etc.) for each atom of interest [29]. Then, analysis of those changes of charge where done.
Additionally, from the NBO outputs corresponding to the intermediates (INT) and transition states (TS) of interest, Wiberg bond indexes (B bond ) [30], were determined and used to calculate the bond evolution (δB bond ) according to Equation (5), which is zero when the bond has not changed, and the unity when the bond has totally evolved. The synchronicity (Sy) [31] of the step of interest was also computed by using Equation (6).
In these equations, i is the i-th step of the mechanism that represents the step of interest, n is the total number of bonds considered to the analysis (bonds that have not significantly changed are neglected), δB bond k is the bond evolution of the k-th bond, and δB bond average is the average of whole δB bond k .

Conclusions
The detailed mechanism of reaction proposed in the present work is in good agreement with the experimentally reported data, given the fact that the polymerization of asparagine is preferred over the diketopiperazine product during pyrolysis at low temperatures (300 • C). The formation of those non-volatile polypeptides justifies the low yields obtained in the LTT's conditions. The high energetic barrier at the rupture process that represents the most-influent rate-controlling process also corresponds well to the necessity to increase the temperature of pyrolysis to obtain the fragmented, five-membered cyclic products (maleimide and succinimide), which are yielded when the temperature exceeds 500 • C. According to the present calculations, the high free energy barrier of the DKP mechanism is associated with the entropic contributions of the dimerization step.
Supplementary Materials: The following are available online at http://www.mdpi.com/2076-3417/9/22/4847/s1. Supplementary Materials include: (1) An Excel document with all calculations to obtain the matrix Mij, which can be used to calculate the degree of rate constants and thermodynamic rate constants of whichever Gibbs free energy profile; (2) the complete matrix Mjk of both mechanisms of pyrolysis of asparagine, Table S1: Matrix Mjk for Mechanism of formation of 2,5-diketopiperazine, Table S2: Matrix Mjk for Mechanism of polymerization-rupture; and (3) the optimized geometries of all reactants, transition states, and intermediates.