Subscriber access provided by CORNELL UNIVERSITY LIBRARY
Article
Model Simulations of the Thermal Dissociation of the TIK(H) Tripeptide. Mechanisms and Kinetic Parameters +
2
Zahra Homayoon, Subha Pratihar, Edward A. Dratz, Ross Snider, Riccardo Spezia, George L Barnes, Veronica Macaluso, Ana Martin Somer, and William Louis Hase J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.6b05884 • Publication Date (Web): 27 Sep 2016 Downloaded from http://pubs.acs.org on September 28, 2016
Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.
The Journal of Physical Chemistry A is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.
Page 1 of 60
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
Model Simulations of the Thermal Dissociation of the TIK(H+)2 Tripeptide. Mechanisms and Kinetic Parameters
Zahra Homayoon,a Subha Pratihar,a Edward Dratz,b Ross Snider,c Riccardo Spezia,d George L. Barnes,e Veronica Macaluso,d Ana Martin Somer,d and William L. Hasea* a
Department of Chemistry and Biochemistry Texas Tech University Lubbock, Texas 79409-1061 USA
b
Department of Chemistry and Biochemistry Montana State University Bozeman, Montana 59717 USA
c
Department of Electrical and Computer Engineering Montana State University Bozeman, Montana 59717 USA
d
Laboratoire Analyse et Modélisation pour la Biologie et l’Environnement Université d’Evry Val d’Essonne UMR 8587 CNRS-CEA-UEVE Bd. F. Mitterrand, 91025 Evry Cedex, France e
Department of Chemistry and Biochemistry Siena College Loudonville, NY 12211, USA
1 ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 2 of 60
Abstract
Direct dynamics simulations, utilizing the RM1 semiempirical electronic structure theory, were performed to study the thermal dissociation of the doubly protonated tripeptide threonine-isoleucine-lysine ion, TIK(H+)2, for temperatures of 1250-2500 K, corresponding to classical energies of 1778–3556 kJ/mol. The number of different fragmentation pathways increases with increase in temperature. At 1250 K there are only 4 fragmentation pathways, with one contributing 85 % of the fragmentation. In contrast, at 2500 K, there are 61 pathways and not one dominates. The same ion is often formed via different pathways and at 2500 K there are only 14 m/z values for the product ions. The backbone and side chain fragmentations occur by concerted reactions, with simultaneous proton transfer and bond rupture, and also homolytic bond ruptures without proton transfer. For each temperature the TIK(H+)2 fragment probability versus time is exponential, in accord with the RRKM and transition state theories. Rate constants versus temperature were determined for two proton transfer and two bond rupture pathways. From Arrhenius plots activation energies Ea and A-factors were determined for these pathways. They are 62-78 kJ/mol and 2-3 x 1012 s-1 for the proton transfer pathways and 153-168 kJ/mol and 2-4 x 1014 s-1 for the bond rupture pathways. For the bond rupture pathways, the product cation radicals undergo significant structural changes during the bond rupture as a result of hydrogen bonding, which lowers their entropies and also their Ea and A parameters relative to those for C-C bond rupture pathways in hydrocarbon molecules. The Ea values determined from the simulation Arrhenius plots are in very good agreement with the reaction barriers for the RM1 method used in the simulations. A preliminary simulation of TIK(H+)2 collision induced dissociation (CID), at a collision energy of 13 eV (1255 kJ/mol), was also performed to compare with the thermal dissociation simulations. Though the energy transferred to TIK(H+)2 in the collisions is substantially less than the energy for the thermal excitations, there is substantial fragmentation as a result of the localized, non-random excitation by the collisions. CID results in different fragmentation pathways with a significant amount of short time nonstatistical fragmentation. Backbone fragmentation is less important and side chain fragmentation more important for the CID simulations as compared to the thermal 2 ACS Paragon Plus Environment
Page 3 of 60
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
simulations. The thermal simulations provide information regarding the long time statistical fragmentation.
*
[email protected], 806-834-3152 3 ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 4 of 60
I. INTRODUCTION Many studies have addressed the chemistry of peptide dissociation, i.e. what are the reaction pathways, activation energies, product ion structures, etc., and representative reviews of this work are in references.1-12 Despite such a large number of studies many fragmentation pathways are not definitively understood and classified.13-18 In addition, during the fragmentation process the peptide sequence can be lost due to scrambling phenomena.19,20 On the basis of this process, there is the possibility of forming intermediate cyclic structures that seem to be the precursors of large macrocycle formation, whose fragmentations give rise to the loss of sequence information.14 Another source of complexity in the fragmentation of a peptide ion comes from the appearance of fragments due to dissociation at two points along the backbone, often including the loss of an internal CO molecule.2,22 There is a need to efficiently identify proteins from mass spectrometry data that contain amino acid sequences not found in current databases. This is significant since methods that strictly rely on sequence databases are prone to flaws in gene finding. Studies have shown that 40-70% of high signal/noise MS/MS spectra cannot be matched to predicted protein spectra (or are misidentified) by widely used search engines.23,24 One of the problems is the post-fragmentation sequence scrambling described above. Another is that peptides do not fragment as expected from simple models (one case is the aforementioned dissociation at two points along the backbone) and are confused with spectra that arise from other peptides. It is likely that this occurs when amino acids undergo side chain fragmentation, along with the usual peptide bond cleavages. This is particularly critical in the case of post-translational modifications.25-27 The conventional theoretical model for MS/MS collision-induced dissociation (CID) spectra of protonated peptide ions, peptide-(H+)n, is that the collision transfers vibrational energy E to the peptide ion and this energy is randomized within the ion via rapid intramolecular vibrational energy redistribution (IVR). RRKM theory28 may then be used to calculate the ion’s dissociation rate constant k(E) and branching between the different fragmentation pathways.29-36 Even if the assumption of rapid IVR is correct, there are uncertainties in applications of the RRKM model. Rotational and vibrational energy have much different effects on the RRKM unimolecular rate constant.28,37 Though 4 ACS Paragon Plus Environment
Page 5 of 60
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
chemical dynamics simulations of CID have shown that a large fraction of the energy transfer to the ion may be to rotation, i.e. as high as ~ 50%,38-40 it is usually assumed the collision energy transfer is to ion vibration.29,36 The importance of energy transfer to rotation is expected to depend upon the structure of the peptide ion as well as the collision partner. In simulations of the collisional excitation of Aln clusters, a large fraction of the energy transfer was to rotation for planar clusters, while transfer to rotation was negligible for spherical clusters.38 A similar result was found for collisional excitation of peptides, where rotational excitation was larger for beta-sheet as compared to folded structures.39 Large rotational energy transfer was also found for protonated urea CID.40 Thus, for some peptide-H+ CID, an assumption that the majority of the collisional energy transfer is to ion vibration may not be valid. This would affect the RRKM analysis and may be critical, since rotation excitation often leads to product formation not predicted by vibration excitation.41 In addition, though anharmonicity may be important,42,43 most RRKM calculations of CID assume harmonic state counting in determining k(E). The anharmonicity may arise from multiple minima, conformers on the potential energy surface (PES)42 and anharmonic molecular vibrations.43 Both of these effects are directly included in a chemical dynamics simulation of unimolecular dissociation (see Section II.B). Furthermore, the CID collisions excite the peptide ions non-randomly and dissociation may occur before complete IVR.44 This has been observed in a number of CID simulations.45-52 Some of this dissociation may be prompt and occur during the CID collision, dynamics referred to as “shattering”.53 Such non-statistical CID dynamics, including shattering, have been characterized in experiments and chemical dynamics simulations for CH3SH+,46,48,54 protonated urea,40,41 and protonated uracil,55,56 and simulations of Cr+(CO)6 and [Li(uracil)]+,47,50 both also studied experimentally.57,58 NonRRKM shattering dynamics have also been observed in simulations53,59-61 and experiments62,63 of surface-induced dissociation (SID). However, shattering was not observed in a simulation of protonated glycine CID.53 A theoretical approach for obtaining MS/MS spectra of CID is to determine the short-time fragmentation pathways from chemical dynamics simulations and for the 5 ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 6 of 60
longer time dynamics, after complete IVR, one may assume the ions dissociate statistically and use the collisional energy transfer distributions obtained from the simulations to determine the vibrational and rotational energy distributions of the undissociated ions.64 The combination of chemical dynamics simulations and statistical unimolecular calculations, to interpret mass spectrometry experiments, was recently discussed for electron ionization mass spectrometry.65 Such chemical dynamics simulations of CID are limited in the time-scale that may be sampled and, thus, generally only fast processes are observed. This limits its use to understand the full fragmentation, in which statistical decomposition is expected at longer times (the importance of statistical fragmentation depends on the experimental conditions). To treat unimolecular decomposition from an atomistic approach and include the use of the statistical RRKM theory, minima and transition states on the potential energy surface must be located. This is doable and relatively easy way for small systems, but it becomes nearly impossible for larger systems like polypeptides. Searching of stationary points is generally done manually, and thus for systems with many atoms and many possible fragmentation pathways this becomes almost impossible. Automatic methods for searching stationary points are promising, but at the present time they find hundreds of such points for small systems.66,67 A possible way to directly obtain the fragmentation products and statistical conditions, with their corresponding rate constants and barriers, is the use of chemical dynamics at different temperatures with random, i.e. statistical, conditions for the unimolecular reactant. This allows performing simulations at high temperatures (at which the system can dissociate in time-scales affordable by simulations) to obtain the temperature dependence of rate constants.43,68,69 In the work presented here chemical dynamics simulations are described for the thermal dissociation of the doubly protonated tripeptide threonine-isoleucine-lysine, TIK(H+)2. The primary structure of TIK(H+)2 is depicted in Figure 1 along with the usual nomenclature used in studying backbone fragmentations.3 The results obtained from the simulations are the fragmentation pathways, their rate constants versus temperature, and their Arrhenius parameters. For a dissociating molecule of the size of TIK(H+)2, the excitation temperature may be associated with a specific energy. Thus, the simulations provide information regarding the mechanisms and rate constants for the statistical 6 ACS Paragon Plus Environment
Page 7 of 60
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
dissociation of TIK(H+)2 versus energy. Furthermore, it will be possible to compare the thermal dissociation and collision dynamics results to determine the importance of nonstatistical fragmentation for the collision activation, and the role of the activation process in the ensuing fragmentation dynamics. There are several reasons TIK(H+)2 was chosen for study. First it is a particularly interesting peptide ion to study because of its side chains and their possible fragmentation dynamics.4,70 Also, in future work the fragmentation dynamics of doubly protonated threonine-leucine-lysine, TLK(H+)2, will also be studied and compared with those for TIK(H+)2. Since isoleucine and leucine have identical masses, but different side chain structures, knowing their fragmentation dynamics will be informative to protein identification when these side chain losses occur. There is limited experimental information regarding side chain fragmentation for I and L, but there is some evidence in the literature4 for different side chain fragmentation for I and L. To interpret biological MS/MS experiments for peptides containing I and L it is important to establish their side chain fragmentation mechanisms. Thus, the likelihood of a particular fragment being lost from I or L needs to be better understood. These fragmentation probabilities for I and L may depend on the neighboring amino acids, so this positional context also needs to be explored. A complete positional characterization will require fragmentation studies of X-I-Y and X-L-Y (where X/Y is a peptide or amino acid), which means there are a total of 800 possible tripeptides that contain isoleucine or leucine in the center. To reduce the number, we first choose as X threonine that has well documented side chain losses.4,70 For Y we have chosen lysine because it can be easily protonated. It is at the N terminal of many tryptic peptides commonly used for protein identification, and has a long side chain that provides for interesting dynamics. The study of a doubly-protonated tripeptide guarantees that only those with a lysine or arginine are likely.
II. COMPUTATIONAL METHODOLOGY A. Electronic Structure Theory and Optimization of Conformers The RM1 semiempirical electronic structure method71 was used for the majority of the direct dynamics simulations reported here. RM1 was chosen for the simulations, 7 ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 8 of 60
since this method is expected to give more accurate barrier heights and reaction enthalpies for organic and biochemical reactions than other semiempirical methods.71 RM1 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.71 RM1 has been used successfully in previous simulations of peptide-H+ dynamics.60,61,72-75 A small number of trajectories were calculated with the AM1 semiempirical method76,77 to compare with the RM1 results. The use of DFT or an ab initio electronic structure method for the direct dynamics simulations is not computationally practical for the system studied here. As shown in previous work,61 RM1 barrier heights for gly-H+ and gly2-H+ dissociation, forming singlet products, are in overall good agreement with values determined with the B3LYP, QCISD(T), MP2, and CR-CC(2,3)78,79 electronic structure theories. AM1 and RM1 energetics are similar for gly-H+ dissociation.61 CR-CC(2,3) gives energies for homolytic bond breakage reactions in excellent agreement with experiment. For homolytic bond ruptures of the gly2-H+ backbone to form doublet radicals, RM1 gives dissociation energies higher that CR-CC(2,3).61 Energy transfer and shattering fragmentation dynamics for gly-H+ + diamond{111} collisions have been compared with AM1 and MP2 direct dynamics,53,80 and these two methods give similar results.
A
list
of
RM1
applications
is
maintained
at
http://www.rm1.sparkle.pro.br/research-topics. To choose initial conditions for the simulations, it was necessary to find the potential energy minimum for TIK(H+)2. A molecular dynamics simulation was performed to find a distribution of structures for TIK(H+)2, which were optimized. The goal was to find the minimum energy conformer, but two additional low energy conformers were found. These optimizations were performed with both RM1 and DFT/B3LYP/6-31G* and both methods give similar results. The structures and energies for the three conformers are given and compared in Figure 2. In further RM1 optimizations an additional conformer was found, with a higher energy of 38.4 kJ/mol and a structure somewhat like that of the highest energy conformer of Figure 2 (see Figure S1 of Supporting Information).
There are important hydrogen bonds for
conformers in Figure 2. For the “blue” hydrogen bond of the minimum energy 8 ACS Paragon Plus Environment
Page 9 of 60
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
conformer, the RM1 and B3LYP bond lengths are 1.68 and 1.79 Å, respectively. For the “blue” and “red” hydrogen bonds of the next in energy conformer, the respective RM1:B3LYP lengths (Å) are 1.70:1.75 and 2.52:1.85. For the “red” hydrogen bond of the highest energy conformer, the respective RM1 and B3LYP bond lengths are 1.74 and 1.92 Å. B. Direct Dynamics Simulations of Thermal Activation The direct dynamics simulations were performed with the VENUS/MOPAC software package, for which the VENUS chemical dynamics computer program81,82 is interfaced with the MOPAC electronic structure code.83 As shown84,85 and applied43,85,86 in previous work, if s ≈ (s-1) for the unimolecular reactant’s s vibrational modes and the unimolecular dissociation energy Eo is much less than the reactant’s energy E, i.e. Eo/E
