Bond dissociation energies of X–H bonds in proteins

Knowledge of reliable X–H bond dissociation energies (X = C, N, O, S) for amino acids in proteins is key for studying the radical chemistry of proteins. X–H bond dissociation energies of model dipeptides were computed using the isodesmic reaction method at the BMK/6-31+G(2df,p) and G4(MP2)-6X levels of theory. The density functional theory values agree well with the composite-level calculations. By this high level of theory, combined with a careful choice of reference compounds and peptide model systems, our work provides a highly valuable data set of bond dissociation energies with unprecedented accuracy and comprehensiveness. It will likely prove useful to predict protein biochemistry involving radicals, e.g., by machine learning.


Introduction
The radical chemistry of proteins is of fundamental importance in biochemistry. Radicals are involved in many enzymatic reactions. Also, oxidative damage to proteins is generally initiated by radicals. For both, knowledge of reliable X-H bond dissociation energies (BDEs, X = C, N, O, S) for amino acids in proteins is important. X-H BDEs allow to predict which hydrogen atom abstractions will be thermodynamically possible, and which lesion sites can be repaired by cellular antioxidants. Preferred sites of radical attack and hydrogen atom abstraction are those associated with the lowest energy barriers, and, according to the Bell-Evans-Polanyi principle, 1,2 these reactions are typically the most exothermic. Therefore, X-H BDEs can also provide a valuable input for machine learning models for the elucidation of the kinetics and mechanisms of radical reactions in proteins.
The BDE is dened as the enthalpy required to homolytically cleave a given bond X-Y: 3

X-H / Xc + Hc BDE = H(Xc) + H(Hc) − H(X-H)
Instead of computing BDEs directly using eqn (1), BDEs can also be calculated using the isodesmic reaction method. 4 An isodesmic reaction is a reaction in which the number of bonds of each bond type is preserved. If an experimental BDE of a reference compound Y-H, BDE(Y-H) exp , is known, one can compute the enthalpy change associated with an isodesmic hydrogen atom transfer (HAT) from a compound of interest X-H to that reference compound Yc, and use it together with BDE(Y-H) exp to obtain the BDE of the bond of interest, BDE(X-H): DH DFT ¼ HðY À HÞ þ HðX Þ À HðX À HÞ À HðY Þ ¼ BDEðX À HÞ À BDEðY À HÞ The isodesmic reaction method is expected to give more accurate BDEs than if they were computed directly because systematic errors of the chosen functional should compensate to some extent if the reference compound is chosen to be electronically similar.
Previously, X-H BDEs for amino acids in proteins have been computed by Moore et al. 5 using the isodesmic reaction method at the B3LYP/6-31G(d) level of theory. This study was an important step forward, as it was the rst detailed catalogue of BDEs for proteins. However, while in the past oen used, B3LYP is one of the less well performing functionals for radical reactions. 6 Also, the set of references used previously 5 lacks proper reference compounds for benzylic C-H, carboxylic O-H, and aromatic N-H bonds. Instead, tert-butyl amine was used as a reference compound for all side chain N-H bonds.
We report X-H BDEs of amino acids in proteins computed using the isodesmic reaction method at the BMK/6-31+G(2df,p) 7 level of theory. We validate our BMK results by comparing the BDEs of smaller amino acids to results obtained at the G4(MP2)-6X 8 level of theory. The BMK functional was chosen because it has been shown to perform well for both, thermodynamics and kinetics of radical reactions. 9-11 G4(MP2)-6X is a composite method that uses BMK/6-31+G(2df,p) geometries and performs almost as well as G4 at lower computational cost. The G4(MP2)-6X energy is obtained as follows: HF/CBS is a complete basis set (CBS) extrapolation of the Hartree-Fock (HF) energy. E corr are frozen-core correlation energies. E corr S-CCSD and E corr S-(T) are scaled CCSD and perturbative triples correlation energies beyond E corr SCS-MP2 /6-31G(d). SCS-MP2 is a modied MP2 method. 12 HLC is a higher-level correction term and E SO is a spin-orbit correction from experiments or accurate calculations. The zero-point vibrational energy (ZPVE) is obtained using scaled BMK/6-31+G(2df,p) frequencies. In a benchmark study of reaction energies for C a -H abstraction in amino acids by HOc, G4(MP2)-6X showed the best performance out of the tested methods. 13 This study builds upon and extends previous results in four aspects: (1) by choosing a level of theory that is very well suited for the treatment of radical reactions, (2) by using a set of reference compounds which correspond better to the local environment in proteins by including more neighbors of the radical center, (3) by using model peptides with C and N termini extended by the next two atoms in the protein backbone, thereby directly modeling the local environment of the amino acid in a protein, and (4) by including BDEs for aspartate, glutamate, dihydroxyphenylalanine (DOPA), hydroxyproline, hydroxylysine, an acetylated N terminus, and an anionic or amidated C terminus, as well as for the lowest-lying tautomer of arginine. The resulting data set, being both reliable and comprehensive, can provide the basis for computationally studying the rich radical chemistry of proteins.

Methods
BDEs in this work were mainly calculated using the hybrid BMK functional 7 with the 6-31+G(2df,p) basis set. For smaller amino acids, BDEs were also computed at the G4(MP2)-6X 8 level of theory. Additionally, calculations were performed using B3LYP/ 6-31G(d) to compare to previous literature results.
Reference compounds and bonds used in the isodesmic reaction method are depicted in Fig. 1. Experimental BDEs were taken from Luo, 14 except for alanine, for which no experimental BDE is available. Instead, the BDE was computed directly using G4(MP2)-6X. Including alanine as a reference compound seemed preferable to the use of another compound because C a radicals are captodatively stabilized, 15 and a reference compound should include this special electronic stabilization if compensation of systematic errors is desired. BDEs of reference bonds computed at the G4(MP2)-6X and BMK/6-31+G(2df,p) levels of theory agree well with experimental values 14 (Table  SI.1 †). The mean absolute errors (MAEs) are 6.2 kJ mol −1 and 8.6 kJ mol −1 , respectively.
In addition to BDEs, for the sake of completeness, bond dissociation free energies (BDFEs) are reported as well. BDFEs were calculated directly, i.e., analogously to eqn (1).
Further experimental details are given in the ESI. † Optimized structures and results of the calculations are available. 26

Results and discussion
The electronic structure of the amino acids and thereby the BDEs are sensitive to their chemical surrounding. Hence, for proteins, adequate capping groups at the N and C termini that model the inuence of the protein backbone must be used.
To assess to what extent the way the protein backbone is modeled inuences the computed BDEs, we compared BDEs of the C a atom with varying capping groups. The way the protein backbone is modelled can be expected to have a particularly strong inuence on the C a -H BDEs, not only because of the close proximity of the C a -H bond to the capping groups, but also because its radical is captodatively stabilized by the backbone amide groups. C a -H BDEs of glycine and alanine were calculated with three different types of capping groups: without, Fig. 1 Reference compounds used to compute BDEs using the isodesmic reaction method. Reference bonds are highlighted and numbered. Triangles indicate values that were computed directly, circles indicate values computed using the isodesmic reaction method. Table 1 BDEs and BDFEs in kJ mol −1 , computed at the BMK/6-31+G(2df,p) level of theory. BDEs were computed using the isodesmic reaction method and the reference bond given in the last column. BDFEs were computed directly. AA, amino acid. with smaller capping groups corresponding to the next backbone atom, and with larger capping groups corresponding to the next two backbone atoms. Calculations were performed using BMK/6-31+G(2df,p) and G4(MP2)-6X, as well as using B3LYP/6-31G(d) for comparison to previous literature results. 5 G4(MP2)-6X BDEs were computed directly, density functional theory (DFT) BDEs were obtained directly and using the isodesmic reaction method. The results are shown in Fig. 2.
Capping alanine with a formyl and an amino group raises the C a -H BDE by around 30 kJ mol −1 compared to the free amino acid by lowering the captodative stabilization by the electron donating and withdrawing groups. This indicates a signicant inuence of the protein backbone on BDEs. Calculations at the G4(MP2)-6X level of theory further show that, although the C a -H BDE of alanine as a free amino acid is lower than that of glycine as a free amino acid, already when considering the next atoms along the peptide backbone the Deviates from a fully extended geometry (f z −164°, j z 91°). e Formyl capping group on N terminus. relative order is reversed. The difference between the two BDEs becomes larger if one includes the next two atoms. This is in agreement with previous calculations at the G3(MP2)-RAD level of theory. 16 The trend might be due to increased steric interactions for the C a alanine dipeptide radical that outcompete the inductive stabilization (a trend that appears reversed for leucine, see Table 1). 17 Therefore, both composite methods come to the conclusion that the C a -H BDE of glycine in a protein should be lower than that of alanine. At the BMK/6-31+G(2df,p) level of theory, this is reproduced for the acetyl and N-methyl capping groups. Hence, all X-H BDEs in this work were computed using these capping groups. The resulting peptide models are shown in Fig. 3. At the B3LYP/6-31G(d) level of theory, the BDE of the alanine dipeptide remains higher than that of the glycine dipeptide even if acetyl and N-methyl capping groups are used, although the exact values are closer to the G4(MP2)-6X results than at the BMK/6-31+G(2df,p) level of theory. Consequently, if one were to use B3LYP, at least the next three atoms along the protein backbone would have to be included. Compared to literature results computed with other composite methods using the isodesmic reaction method, G4(MP2)-6X seems to underestimate C a -H BDEs, but BMK C a -H BDEs computed using the isodesmic reaction method agree very well with composite level literature results (Table SI.2 †). BMK BDEs are also closer to high-level composite methods than DSD-PBE-P86/aug ′ -cc-pVTZ+d BDEs reported by Chan et al. 13 DFT, especially B3LYP, underestimates BDEs compared to G4(MP2)-6X if they are computed directly. Using the isodesmic reaction method, a large part of this deviation vanishes. This highlights the error-canceling power of the isodesmic reaction method.
We note that we were not able to reproduce the C a -H BDEs of glycine and alanine dipeptides reported by Moore et al., 5 or the BDEs of their reference compounds, although our calculations agree well with previous literature results. [18][19][20] Their reported BDEs also differ from previous literature results at the same level of theory and using the same reference compounds by as much as 9 kJ mol −1 , 19 and by over 30 kJ mol −1 compared to a study that uses glycine as a reference compound instead of alanine. 20 Importantly, we and others, 16 using composite methods or the BMK functional and the isodesmic reaction method, nd the C a -H BDE of the glycine dipeptide to be lower than that of the alanine dipeptide, in contrast to Moore et al., who also report overall much lower BDEs for these two radicals.
BMK BDEs are given in Table 1. G4(MP2)-6X BDEs of smaller amino acids can be found in Table 2. In general, BDEs at the BMK/6-31+G(2df,p) level of theory agree well with BDEs at G4(MP2)-6X level of theory (coefficient of determination of the linear t, R 2 = 0.997, mean absolute deviation, MAD = 4.7 kJ mol −1 , Fig. 4). Some X-H BDEs were also computed for model peptides with formyl and amino capping groups. The BDEs for model peptides with larger and smaller capping groups correlate well with R 2 = 0.981 (Fig. SI.1 †), and just the C a -H BDEs deviate unsystematically from the linear correlation, as also discussed above (see also Fig. 2). It is therefore sufficient to use C a -H BDEs to determine adequate capping groups.  For the BDEs given in Table 1, peptides were modeled in a fully extended trans geometry (u z 180°), except for proline and hydroxyproline, where both cis and trans conformers were modeled by setting u, f, j, and u ′ going from the N to the C terminus to approximately 0°, −75°, 160°, 180°, and 180°, −75°, 150°, 180°, respectively. 21 Where applicable, amino acids were considered in their charged and neutral forms. For histidine, the neutral N 1 -H and N 3 -H tautomers were considered. In the gas phase, the N 3 -H tautomer is more stable by 4.9 kJ mol −1 . For arginine, in agreement with previous investigations, 22 out of the ve possible neutral tautomers, the one in which the two terminal N atoms are doubly protonated was found to be most stable, and so this tautomer was considered.
In some cases, initial geometry optimizations resulted in hydrogen bonding between the amino acid side chain and backbone to different extent in radical and non-radical structures. This intramolecular hydrogen bonding is well known, and has also been found in a comprehensive study of the conformational landscape of amino acid dipeptides. 23 For the computation of BDEs, this would lead to energetic artifacts. Intramolecular hydrogen bonding is unlikely to occur in proteins, since the backbone amide groups are involved in secondary structure formation. In aqueous solution, polar side chain groups will also be hydrogen-bonded to water molecules or coordinated to metal ions. Therefore, we aimed at minimizing intramolecular hydrogen bonding and changes in hydrogen bonding pattern. To this end, the same initial structure was used for geometry optimizations of both dipeptides and dipeptide radicals to reect the potential structural constraints imposed by the protein that would make large conformational changes unlikely. Further efforts to this end are described in the ESI. † Yet, in some cases, signicant changes in hydrogen bonding strength or pattern were observed, as speci-ed in further detail in Table 1.  BMK and G4(MP2)-6X BDEs overall agree well with other literature results for the same model peptides (Table SI.3 †). Deviations from the BDEs reported by Zipse and colleagues 17,24 can be attributed to the fact that they use Boltzmann averages. The resulting difference is small for the BDEs in tyrosine, proline, and phenylalanine, since these amino acids have less conformationally exible side chains, but becomes apparent for the BDEs in cysteine, where the side chain is more exible. Deviations from the BDEs reported by Chan et al. 13 for amino acids with polar side chains come from the fact that they have used the lowest energy conformations of the peptides, which involve hydrogen bonding for those amino acids. In this work, hydrogen bonding in general and in particular changes in hydrogen bonding pattern between peptides and peptide radicals were avoided, and by design of our protocol the resulting conformers do not necessarily represent the most stable conformer.
Having a comprehensive data set of BDEs at hand, we analyzed trends in the BDEs across the chemical space covered by the amino acids. We grouped the X-H BDEs according to the reference compound used in the isodesmic reaction and show the distribution of the ten largest groups as well as an example radical for each group in Fig. 5. Since the reference compounds were chosen to correspond to the local chemical environment of the radical, Fig. 5 shows the inuence of the local chemical environment on the BDEs. The observed stability trends agree with known chemical principles: captodative and resonance stabilization leads to low BDEs, as can be seen for C a -H, as well as benzylic C-H and O-H bonds. Secondary C-H bonds and C-H bonds adjacent to a heteroatom have medium high BDEs. Primary C-H BDEs are higher. N-H, O-H, and aromatic C-H bonds are associated with the highest BDEs.
Finally, we sought to evaluate the utility of our data set in the context of machine learning. To this end, we compared our data set to predictions made by ALFABET, a graph neural network for the prediction of BDEs trained on small organic molecules consisting of up to ten heavy atoms. 25 Considering that most dipeptides in our data set consist of more than 10 heavy atoms, ALFABET predictions agree reasonably well with BMK BDEs (R 2 = 0.918, MAD = 11.3 kJ mol −1 , Fig. SI.5 †). There is, however, considerable scatter for lower BDEs, suggesting that ALFABET does not fully capture the chemistry of amino acid radicals. Our data set could therefore prove useful in a transfer learning setting for the radical chemistry of proteins.

Conclusion
BDEs of all X-H bonds in amino acids were computed using the isodesmic reaction method at the BMK/6-31+G(2df,p) and G4(MP2)-6X levels of theory. The BDE values agree well with high-level calculations, and improve upon previous literature results by using a more appropriate level of theory, reference compounds, and peptide model systems.

Conflicts of interest
There are no conicts to declare.