Direct Dynamics Simulations of the Thermal Fragmentation of a Protonated Peptide Containing Arginine.

Arginine has significant effects on fragmentation patterns of the protonated peptide due to its high basicity guanidine tail. In this article, thermal dissociation of the singly protonated glycine-arginine dipeptide (GR-H+) was investigated by performing direct dynamics simulations at different vibrational temperatures of 2000-3500 K. Fourteen principal fragmentation mechanisms containing side-chain and backbone fragmentation were found and discussed in detail. The mechanism involving partial or complete loss of a guanidino group dominates side-chain fragmentation, while backbone fragmentation mainly involves the three cleavage sites of a1-x1+, a2+-x0, and b1-y1+. Fragmentation patterns for primary dissociation have been compared with experimental results, and the peak that was not identified by the experiment has been assigned by our simulation. Kinetic parameters for GR-H+ unimolecular dissociation may be determined by direct dynamics simulations, which are helpful in exploring the complex biomolecules.

Introduction

Research on the protein structure and function is of great significance since proteins play a vital role in all aspects of life activities.[1,2] Mass spectrometry has become a key technology for determination of amino acid sequence information in proteins.[3−5] However, due to the complex composition of proteins, many mass spectrometry data are often not fully understood and classified.[6−8] To identify the amino acid sequence of proteins, it is important to enhance the understanding of protein dissociation patterns. Computational simulations provide a reliable way to study the dissociation dynamics of peptide ions and their fragmentation patterns, providing detailed knowledge of protein composition at the atomic level.[9−11]

Arginine is important in mass spectrometric examinations of proteins since it has a high proton affinity side chain that can effectively capture the mobile proton, resulting in fewer backbone fragmentations for singly protonated peptides containing arginine.[12] Moreover, peptide ions are prone to undergo cyclization and rearrangement reactions, yielding a large amount of fragment ions, resulting in spectra that are complicated and difficult to resolve.[13,14] The presence of arginine may greatly affect fragmentation patterns of peptide ions. Glish et al. found that the cleavage of the peptide bond between arginine and the adjacent amino acid requires more energy than the peptides without arginine and suggested that secondary interactions between the arginine side chain and an adjacent amino acid could be responsible for the effects of arginine on dissociation patterns of peptides.[15] Paizs et al. investigated proton-driven cleavage of the amide bond in arginine-containing peptide ions (GR, n = 2–4) and proposed three new mechanisms, including salt-bridge, anhydride, and imine enol intermediates.[16]

The protonated glycine-arginine dipeptide, abbreviated as GR-H+, is one of the simplest dipeptide ions containing arginine, and investigations of its dissociation and cleavage sites would be very helpful for better understanding of more complicated peptides containing arginine. O’Hair et al.[17] demonstrated that collisional activation of protonated Gly-Arg-H+ and Arg-Gly-H+ gave the same collision induced dissociation (CID) mass spectra by using a quadrupole ion trap instrument. Forbes et al.[14] continued the exploration of the effect of different activation methods on the fragmentation of above two protonated dipeptides and suggested that the variety of experimental conditions may affect the experimental mass spectra. For some conditions, rearrangement of the dipeptide ions to a common structure may occur before fragmentation, leading to deviations in the mass spectra.

Since fragmentation patterns of peptide ions may be influenced by the activation methods, it is important to determine the intrinsic fragmentation mechanisms for peptide ions containing arginine.[18−22] Electronic structure calculations can identify important dissociation pathways and their transition states (TSs).[23−26] However, for large molecules such as peptides, the important dissociation pathways may be difficult, or even impossible, to determine owing to the numerous fragmentation possibilities. Direct dynamics simulations[27,28] provide an efficient approach to interpret the mass spectrometry studies of peptide ion fragmentation. For the simulations, the peptide ions may be excited with thermal vibrational energies randomly distributed among the vibration mode and sufficiently high to allow full fragmentation within a practical simulation time scale. Hase et al. have successfully performed direct dynamics simulations of the thermal dissociation of the doubly protonated tripeptide threonine-isoleucine-lysine ion, TIK(H+)2.[29] Important dissociation pathways of the peptide ion, and their corresponding reaction probabilities, were determined versus temperature. For each temperature, the TIK(H+)2 fragmentation probability versus time was exponential, consistent with Rice–Ramsperger–Kassel–Marcus (RRKM) theory and transition state theory (TST).

The current study primarily focuses on fragmentation patterns of arginine-containing peptide ions GR-H+ by performing direct dynamics simulations. The primary structure of GR-H+, and its corresponding fragmentation nomenclature proposed by Roepstorff and Fohlman,[30] is depicted in Figure . Complete identification of the peptide ions’ dissociation dynamics, including atomistic dissociation mechanisms, fragmentation pathways, and reaction probabilities, was determined by varying the vibrational temperature for the simulations. The competition of side-chain and backbone dissociations is observed and discussed. By determining the unimolecular rate constant k(T) and reaction probabilities for each temperature T, the Arrhenius A-factor (A) and activation energy (Ea) were determined for different fragmentation pathways. Through comparisons of the simulation fragmentation patterns with available experimental data, we expect to establish how the presence of arginine affects the fragmentation dynamics of arginine-containing peptides. This may provide some explanations and insights into the MS spectra of arginine-containing peptides and further enrich the knowledge of unimolecular dissociation dynamics.

Figure 1

Primary structure of GR-H+ and fragmentation nomenclature proposed by Roepstorff and Fohlman.[30]

Results and Discussion

Structure of GR-H+

The GR-H+ potential energy minimum is needed for selecting initial conditions for the direct dynamics trajectories. As shown in Figure , three possible protonated sites, that is, imidic nitrogen of the guanidino group (a), amino of the N-terminus (b), and carbonyl oxygen (c), were considered for performing the optimizations. It was found that configurations with the protonated guanidino group (a) give relatively lower energies than for the other two conformers (b and c), which have energies 167.7 and 189.9 kJ/mol higher, a finding consistent with Paizs’s work for protonated ariginine.[12] By referencing previous works,[17,31] the three lowest-energy GR-H+ conformers for protonation at guanidino (a) were located and are displayed in Figure a. We found that the energy of the conformer II given in ref (31) as the lowest-energy conformer is close to that of conformer III with the difference of 0.4 kJ/mol with the RM1 method, and both of them are higher in energy than conformer I. As shown in Figure a, protonated guanidine of the lowest-energy conformer I has an extended conjugated π bond, which is coplanar with a dihedral angle of 179.1°. The three N–C bond lengths and N–C–N angles are almost the same with the values of ∼1.38 Å and 120.0°, respectively. This is the most stable conformation, with the protonation of the guanidino group chosen as the initial configuration for the direct dynamics simulations. In previous works for both TIK(H+)2 and TLK(H+)2[32,33] collisional simulations, Hase and co-workers found that the lowest-energy conformers and the next lowest-energy conformers both gave statistically the same fragmentation dynamics. The reason is because the small difference in the potential energy minima of the two conformers is overwhelmed by their zero-point energies (ZPEs), which are included in the trajectory simulations. For conformers I–III in Figure in our work, the harmonic ZPEs are 183.1, 183.0, and 183.2 kcal/mol, respectively, which are much larger than the energy differences between them.

Figure 2

Conformers and their relative energies for GR-H+ obtained by RM1: (a) three lowest-energy conformers for GR-H+ with protonation at the imidic nitrogen of the guanidino group, (b) conformer for protonation at the amino of N-terminus, and (c) conformer for protonation at the carbonyl oxygen. Energies are given in kJ/mol and distance in angstrom.

Fragmentation Probability and Total Rate Constant

There are 172 dissociating trajectories out of the 400 total trajectories, and they may be classified by 66 different primary fragmentation pathways. Here, we focus on the primary dissociation pathways and do not consider secondary dissociations of primary dissociation products. This latter topic would be of interest for future simulations. The number of primary dissociation pathways tends to gradually increase as the temperature is increased, with values of 16, 20, 22, and 40 at 2000, 2500, 3000, and 3500 K, respectively. Correspondingly, the total fragmentation probability increases with the increase in temperature with values of 0.20, 0.31, 0.47, and 0.74, as shown in Figure . Here, we define pathways, with a fragmentation probability greater than 3%, as principal dissociations. There are 14 principal dissociation pathways and 52 minor pathways. The reaction probabilities for the principal decomposition pathways grow linearly with the increase in temperature, as presented in Figure .

Figure 3

Plot of reaction probability Pr(b) relative to the total trajectories for the unimolecular dissociation reaction of GR-H+ versus temperature (T). The black solid line represents the total reaction probability. The red and green dashed lines are for principal dissociation pathways and minor pathways, respectively. The blue and orange dotted lines are for backbone fragmentation and side-chain fragmentation, respectively.

To discuss the fragmentation patterns of GR-H+ in detail, it is necessary to analyze both backbone and/or the side-chain fragmentations. The probabilities of principal paths involving side-chain fragmentation are 0.03, 0.11, 0.20 and 0.28 for the temperatures of 2000, 2500, 3000, and 3500 K, respectively, as shown in Figure . For backbone fragmentation, the probability increases from 0.09 to 0.19 with T increasing from 2000 to 3500 K. The importance of both backbone and side-chain fragmentation grows with the increase in temperature.

The total unimolecular rate constant at each temperature may be determined by fitting N(t)/N(0) versus time through the equation N(t)/N(0) = exp(−kt), where N(t)/N(0) is the fraction of reactants surviving at time t.[34] The slope of ln[N(t)/N(0)] versus t in Figure gives the unimolecular rate constant k(T) values, which are (6.00 ± 0.27) × 1010, (2.42 ± 0.12) × 1011, (8.73 ± 0.30) × 1011, and (1.65 ± 0.04) × 1012 s–1 for T of 2000, 2500, 3000, and 3500 K, respectively. The linear fit of ln k(T) versus 1/T as presented in Figure yields the Arrhenius parameters, which are Ea (activation energy) = 131.1 ± 6.3 kJ/mol and A (pre-exponential factor) = (1.59 ± 0.46) × 1014 s–1.

Figure 4

Plot of ln[N(t)/N(0)] versus time at different temperatures for the unimolecular dissociation reaction of GR-H+.

Figure 5

Natural logarithm of the overall rate constant (s–1) for GR-H+ dissociation plotted versus 1/T (1 × 10–4 K–1). The Arrhenius parameters are A = (1.59 ± 0.46) × 1014 s–1 and Ea = 131.1 ± 6.3 kJ/mol. The R value for the fit is −0.993. The blue left solid diamonds are the values obtained from the rate constant for GR-H+. The red straight line corresponds to the fit using the Arrhenius equation.

Principal Dissociation Channels and Their Barriers

To focus on the decomposition mechanisms, the number of trajectories for each principal path relative to the total number of dissociating trajectories is depicted in Figure as a function of temperature, and the corresponding proportions of dissociation pathways are given in Table S1 of the Supporting Information. The 14 principal paths and their corresponding transition states are described in Figures and 8. As shown in Figure , the probability of most of the paths displays an ascending trend with increasing temperature. At the higher temperatures of 3000 and 3500 K, path 4 is dominant, while paths 6 and 1 are dominant at 2000 and 2500 K, respectively. Paths 1–5 in Figure involve side-chain fragmentation, and they account for 36.0 ± 3.7% of the total dissociation trajectories. The proportion of backbone breakage is 33.7 ± 3.6%, suggesting that the former mechanism is even comparable to the latter. This phenomenon is consistent with the mobile proton model, in which sequestration of proton by the highly basic side chain leads to the less efficient backbone fragmentation. This result is also consistent with the experimental observation, as discussed below.

Figure 6

Proportions of 14 principal dissociation pathways versus temperatures of 2000, 2500, 3000, and 3500 K. The proportions of each pathway to total dissociation trajectories are listed in the picture. The column length of different colors represents the reaction possibility of each temperature, which is either 3500 K(dark blue), 3000 K(blue), 2500 K(light blue), or 2000 K(gray).

Figure 7

Mechanisms for dissociation pathways 1–5 together with their corresponding transition states: (a) reactions producing NH3, (b) reactions producing NH=C+—NH2, and (c) reactions producing neutral guanidine and protonated guanidine. The barrier heights for the transition states are given by RM1 in kJ/mol. The arrows correspond to the direction of proton transfer.

Figure 8

Mechanisms for dissociation pathways 6–14: (a) reactions producing a1 + x1+, (b) reactions producing a2+ + x0, and (c) the three-body decomposition mechanism. The arrows correspond to the direction of proton transfer.

Side-Chain Fragmentation

Reactions Producing NH3 (Pathways 1 and 2)

As shown in Figure a, both paths 1 and 2 are side-chain fragmentation mechanisms. Proton transfer on the guanidino group leads to loss of NH3, leaving either the —NH+=C=NH or —N=C=NH2+ moiety behind. Path 1 is the most preferred event for RG-H+ dissociation, which accounts for 12.2 ± 2.5% of the total dissociation, in contrast to 3.5 ± 1.4% for path 2, as depicted in Figure . In addition, path 1 occurred at each temperature, whereas path 2 only occurred at the higher temperatures of 3000 and 3500 K. The transition states for these two pathways were calculated using the RM1 method, and their optimized structures are given in Figure a. The TSs for both paths 1 and 2 represent a concerted process, with hydrogen transfer and C–N2 bond breakage occurring simultaneously. For the TSs, the N–H and C–N2 bonds elongated to ∼1.41 and 1.49 Å, respectively, compared to the corresponding distances of 1.01 and 1.36 Å for the minimum of GR-H+. The calculated potential energy barriers for these two pathways are 176.9 and 173.6 kJ/mol.

Reactions Producing NH=C+—NH2 (Pathway 3)

As shown in Figure b, instead of C–N2 bond cleavage as for pathways 1 and 2, the C–N1 bond breaks in this mechanism. Together with hydrogen transfer from N3 to N1 atom, the product NH=C+—NH2 cation is generated, which has the smallest m/z value of 43 among all major products. In addition, it is worth noting that this pathway occurred at all temperatures except for the lowest 2000 K. The TS structure for path 3 has a barrier of 179.8 kJ/mol and exhibits similar characteristics as those for paths 1 and 2 (Figure b).

Reactions Producing Neutral Guanidine and Protonated Guanidine (Pathways 4 and 5)

As described in Figure c, proton transfer is not involved in path 4, and the arginine residue loses the neutral guanidine directly, generating a cation product with the m/z value of 173. By contrast, for pathway 5, the hydrogen atom on the side chain migrates to the N atom in the guanidino group, forming protonated guanidine with the m/z value of 60. As shown in Figure , the contribution of path 4 is 9.3 ± 2.2%, which is larger than that of path 5 (5.2 ± 1.7%), which is probably due to the fact that guanidine is more stable if not protonated. Moreover, path 4 only appears at a higher temperature, and it is the most preferred mechanism at the temperatures of 3000 and 3500 K qualitatively.

Backbone Fragmentation

Reactions Producing a1 + x1+ (Pathways 6–9)

Paths 6–9 belong to the a1-x1+ fragmentation type, and all are concerted reactions with C–C bond cleavage on the backbone of glycine’s N-terminus accompanied by proton transfer. The difference of these four mechanisms is that the hydrogen comes from different groups, either the backbone or the guanidino, leading to two similar products with m/z values of 202 and 203, as presented in Figure a. Path 6 has 8.7 ± 2.2% of the dissociation, and in contrast, the respective contributions of paths 7–9 are 2.9 ± 1.3, 2.3 ± 1.1, and 2.3 ± 1.1%, respectively, much lower than that for path 6, as shown in Figure . This suggests that, for backbone fragmentation, the hydrogen on the backbone migrates more efficiently than that on the side chain since the former proton is closer to the fragmentation site.

Reactions Producing a2+ + x0 (Pathways 10–13)

As shown in Figure b, the products of paths 10–13 have the same m/z value of 186, although the proton comes from different groups. The C–C breakage occurs on the carboxyl of the arginine C-terminus, leaving one HCOOH neutral molecule and the HO–C̈–OH diradical, and an ion with the same m/z value for all four pathways. Formation of HCOOH derives from proton migration on a C atom of the backbone, while HO–C̈–OH is formed by hydrogen transfer from a N atom. Pathway 10 has the highest reaction probability when compared with pathways 11–13, as shown in Figure , again testifying to the more efficient mobility of a backbone proton.

Three-Body Decomposition Mechanism (Pathway 14)

Figure c depicts path 14, which is a three-body decomposition mechanism and contains amide bond cleavage with the probability of 4.7%. A proton on guanidino migrates and attacks the amide N, resulting in the simultaneous occurrence of a1-x1+- and b1-y1+-type backbone fractures, releasing CO and two other products.

Comparison with Experimental Results

The simulation m/z mass spectra for thermal GR-H+ decomposition at 2000, 2500, 3000, and 3500 K are shown in Figure . All m/z values are identified, and the height of the peaks characterizes the reaction probabilities. The total number of m/z values for both principal and minor pathways is 8, 10, 14, and 21 at 2000, 2500, 3000, and 3500 K, respectively. These values do not necessarily correspond to the number of pathways since different pathways can produce the same ion and one pathway can yield multiple ions. Consequently, even though the total number of dissociation patterns is 57, the number of total m/z values is 24.

Figure 9

Theoretical primary m/z mass spectra of GR-H+ trajectories for thermal dissociation at 2000, 2500, 3000, and 3500 K, respectively.

Forbes et al.[14] and O’Hair and Farrugia[17] have experimentally investigated dissociation of protonated GR-H+ by carrying out cone-voltage CID, quadrupole TOF CID, and quadrupole ion trap CID. Mass spectra from the three CIDs show that the main m/z peaks are 215, 214, 175, 173, 158, 100, and 70. Among them, m/z 215 corresponds to dissociation paths 1 and 2, which are the dominant pathways in both the experiments and simulations, resulting from loss of neutral NH3. These reaction pathways are also found in Bythell et al.’s[13] theoretical calculations for fragmentation of singly protonated peptides with N-terminal arginine. The ion m/z 214 in the simulations corresponds to a pathway releasing H2O. It is not a major pathway but is present at both 2500 and 3500 K. It is speculated in the experiment that m/z 175 refers to amide bond cleavage yielding b1-y1, which is consistent with the simulation result. The m/z values of 173 and 60 are related to pathways 4 and 5, which give rise to the loss of neutral and protonated guanidine from the side chain. From the experiments, the m/z 158 peak is thought to arise from the loss of NH3, but its structural formula could not be determined. The simulations suggest that it is OH—CO—CH=CH—(CH2)2—HN—C+—(NH2)2, formed by c1-z1 cleavage. The m/z 203 (x1+) is abundant in the simulations but not observed in experiments of O’Hair and Farrugia.[17] This is because the ion of m/z 203 is metastable and tends to dissociate with more simulation time. As an example, the secondary dissociations of pathway 6, the x+ ion dissociation to form ions of m/z = 175 and 60, are given in Figure S1 of the Supporting Information, in which both ions are observed in the experiments.

Overall, the m/z values obtained from the simulations and experiments are consistent, but there are still some experimental ion signals not found in the theoretical mass spectra. Two reasons may be responsible for the inconsistency. One is that the primary dissociations were considered in current simulations and some ions observed in the experiments may result from the secondary dissociations. The other reason is that more complex pathways may occur and new ions could be observed when the trajectories run enough long time. Our simulations run within several picoseconds, in contrast to the microsecond timescale in the experiments, which may result in the disappearance of some ions in theory.

Conclusions

Different unimolecular dissociation pathways and their probabilities for the protonated glycine-arginine dipeptide were investigated by direct dynamics simulations at different vibrational temperatures of 2000, 2500, 3000, and 3500 K. The calculated thermal dissociation results have been compared with CID experiments, and the corresponding m/z signals in the experiment have been assigned. Some key conclusions are listed as follows:

Fourteen principal dissociation pathways were found and classified into side-chain and backbone fragmentations, and their reaction possibilities tend to ascend with increasing temperature. Consisting with the mobile proton model, the backbone fragmentation is less favored for the singly protonated peptides containing arginine due to the location effect on the proton of guanidine group.

The side-chain breakage is mainly the partial or complete loss of the guanidino group, and the departure of NH3 is a dominant pathway in both the experiment and simulation corresponding to the m/z 215 peak. The backbone fragmentation mainly consists of the cleavage at positions a1-x1+, a2+-x0, and b1-y1+. The dominant a1-x1+-type fragmentation is the concerted reaction with C–C bond cleavage and proton transfer occurring simultaneously.

The rate constants and Arrhenius parameters can be determined by direct dynamics simulations, which is meaningful in exploring the larger reaction systems.

Computational Methodology

Electronic Structure Theory

Since the gradient and potential energy for the trajectories are directly obtained from electronic structure theory, an appropriate method needs to be chosen for the simulations. It is common to use DFT or MP2 for small molecules, while for larger molecules such as peptides, it is more practical to use semi-empirical Hamiltonians, such as AM1, PM3, or RM1, which have been widely applied to macromolecular systems.[35−40] In this work, the RM1 method[36] was chosen for the direct dynamics simulations because it was parameterized with a training set of 1736 molecules, including both neutral and protonated amino acids and representative examples of dipeptides with α-helix and β-sheet conformations,[41] and thus, it is expected to provide more accurate barrier heights and reaction enthalpies for organic and biochemical reactions.[29,35−38] The transition states (TSs) are located for the primary dissociation pathways 1–3 in the light of the trajectory results. It is noteworthy that our simulations show that most of the dissociation processes are complex, and the secondary or tertiary dissociations often occur. In this paper, since we only focus on the primary dissociation pathways, the TSs usually correspond to a simple direct process of proton transfer. In a previous work on the surface-induced dissociation (SID) of protonated glycine, gly-H+, with a diamond surface, it was found that AM1 and MP2/6-31G(d) yield similar SID dynamics.[9,42]

Direct Dynamics Simulations of Thermal Activation

Direct dynamics simulations[27] were performed with the VENUS/MOPAC[43−45] software package. As discussed and applied previously,[46−49] for a molecule with s harmonic oscillators, the classical RRKM and TST rate constants become identical when the unimolecular dissociation energy E0 is much less than the reactant’s energy E and s is large (i.e., E0/E ≪ 1 and s ≈ s – 1), resulting in the relationship of E = skBT. We refer our computational results to experimental results reported by Forbes et al., where CID energies were 12, 20, 25, and 40 eV, respectively.[14] To compare with the experiment, here, simulations of the thermal dissociation of GR-H+ were carried out at temperatures of 2000, 2500, 3000, and 3500 K. For each temperature, the initial vibrational energy was randomly distributed among the vibration mode, where s for the GR-H+ dipeptide is 96. Therefore, the total vibrational energy of the ion was around 17, 21, 25, and 29 eV for 2000, 2500, 3000, and 3500 K, respectively. These energies are in between the experimental collision energies of 12 and 40 eV.[14] It should be noted that less than 100% of the experimental energies can be converted into internal energy of the ion, so the simulation results can be compared with the experimental observations qualitatively. We expect to provide an understanding of mechanisms for the dissociations of arginine-containing peptides in this work. A 300 K rotational energy determined by RT/2 was added to each principal axis of rotation for the peptide ion, which was randomly rotated about its Euler angles.

Direct dynamics simulations were performed using RM1. Different integration times of 27, 7, 3.3, and 1.6 ps were used for the temperature of 2000, 2500, 3000, and 3500 K, respectively. Since the lifetime of unimolecular reaction decreases with the increase in temperature, the required integration time becomes shorter with the increase in temperature. For each temperature, 100 trajectories were calculated. Hamilton’s equations of motion were numerically integrated with a sixth-order symplectic algorithm.[50,51] Different integration step sizes were examined to obtain an acceptable energy conservation, and step sizes of 0.1, 0.05, 0.01, and 0.01 fs were used for the above respective temperatures.