Hydrogen bond donors and acceptors are generally depolarized in α‐helices as revealed by a molecular tailoring approach

Hydrogen‐bond (H‐bond) interaction energies in α‐helices of short alanine peptides were systematically examined by precise density functional theory calculations, followed by a molecular tailoring approach. The contribution of each H‐bond interaction in α‐helices was estimated in detail from the entire conformation energies, and the results were compared with those in the minimal H‐bond models, in which only H‐bond donors and acceptors exist with the capping methyl groups. The former interaction energies were always significantly weaker than the latter energies, when the same geometries of the H‐bond donors and acceptors were applied. The chemical origin of this phenomenon was investigated by analyzing the differences among the electronic structures of the local peptide backbones of the α‐helices and those of the minimal H‐bond models. Consequently, we found that the reduced H‐bond energy originated from the depolarizations of both the H‐bond donor and acceptor groups, due to the repulsive interactions with the neighboring polar peptide groups in the α‐helix backbone. The classical force fields provide similar H‐bond energies to those in the minimal H‐bond models, which ignore the current depolarization effect, and thus they overestimate the actual H‐bond energies in α‐helices. © 2019 The Authors. Journal of Computational Chemistry published by Wiley Periodicals, Inc.


Introduction
The hydrogen bond (H-bond) is one of the major factors that build the macromolecular structures of proteins, nucleic acids, and their complexes. In particular, pair-wise H-bonds in protein backbones are essential to form their characteristic threedimensional (3D) structures based on their ordered secondary structures, α-helices, and β-sheets. Therefore, their structural energies should be correctly computed for analyses and predictions of protein 3D structures.
The individual force fields used in classical molecular dynamics (MD) simulations show particular preferences to produce α-helices and β-sheets. [1][2][3] This phenomenon is usually not a problem for simulations of rigid globular protein structures, but it becomes a critical issue for folding simulations of flexible disordered regions [4,5] and long loops between secondary structures, [6,7] to understand the functionally important conformational changes that occur as allosteric effects or induced folding upon ligand binding. [4][5][6][7] Many attempts have been made to overcome this problem, by improving or rearranging the torsion energies, [8][9][10][11] and by developing polarized charge models. [12,13] However, these preferences have remained unclear, because their actual origins are unknown.
In our previous study, [14] we computed the conformation energies of the secondary structures formed by peptide fragments, using several quantum chemical (QM) methods: the Hartree-Fock (HF) method, the second-order Møller-Plesset perturbation theory (MP2), and the density functional theory (DFT). Consequently, we found that a high-quality DFT method including van der Waals interactions, B97D/6-31+G(d), was comparable to the MP2 method, which is reliable but time-consuming, for the Ace-(Ala) n -Nme system in vacuo. [14] Using this DFT method, the energies of parallel and anti-parallel β-sheets can be approximated more or less by the classical force fields, AMBER ff99SB, [15] but those of the α-helical structures are significantly different. The differences were suggested to originate from the electrostatic energies associated with the H-bonds. [14] In this article, by using the molecular tailoring approach (MTA) [16] with the DFT B97D/6-31+G(d) method, we dissected the individual interaction energy associated with each H-bond, to form typical α-helix backbones with different lengths. To analyze the origin of the H-bond interaction energy in an α-helix, we designed additional simplified models: a minimal Hbond (MH) model, composed of only the atoms forming a single H-bond, and a single-turn (ST) model, composed of three successive alanine residues, (Ala) 3 , in the α-helix capped by acetyl and N-methyl groups at the N-and C-termini, respectively. For those models, the individual H-bond energies were also computed by using MTA with the same DFT method, and the differences in the H-bond energies and the electronic structures among the complete α-helix and several models were analyzed. Finally, we discuss the putative reason underlying the secondary structure preferences in the classical force fields used in molecular mechanical (MM) calculations.

Materials and Methods
α-Helical structure, single turn, and minimal H-bond models The α-helix models were first constructed by using 3-to 8-mer poly-alanine amino acids capped with an acetyl group (Ace) and an N-methyl amide group (Nme), denoted as Ace-(Ala) n -Nme, with the uniform (φ, ψ, ω) backbone angles for each residue: φ = −57 , ψ = −47 , and ω = 180 . Here, n is from 3 to 8. Each structure was optimized in vacuo by energy minimization of the electronic state, while maintaining the backbone angles as mentioned below.
There are one to (n-2) backbone hydrogen-bonds (H-bonds) in the optimized Ace-(Ala) n -Nme (3 ≤ n ≤ 8) α-helical structures, which are denoted here as "AH models," between the backbone carbonyl group (─C O) and the amide group (─NH) from the N-to C-terminus. They are denoted as AHn-1 to AHn-(n-2), respectively (Fig. 1A), and their H-bond energies were individually computed by using MTA with the DFT method.
To analyze the origin of the H-bond interaction energy in the α-helix, we designed a minimal H-bond model (MH model), which is composed of two separated N-methylacetamide molecules, mimicking a single H-bond between the i-th and (i+4)-th residues (Fig. 1B). In addition, we designed a ST model, which is composed of three successive alanine residues, (Ala) 3 , in the α-helix capped by acetyl and N-methyl groups at the N-and Ctermini, respectively (Fig. 1C). The atom positions in the MH and ST models were the same as those in the corresponding AH models, except for the N-and C-terminal capping groups. For two models, the individual H-bond energies were computed by using MTA with the DFT method, in the same manner as for the AH models, and the energy differences among the H-bonds in those models were analyzed.

Theoretical calculations
All calculations including the structure optimization were performed on the above atomic models with the Gaussian09 program package [17] with the DFT B97D/6-31+G(d) method, which can correctly reflect the van der Waals interactions, and this method was confirmed to be comparable to the MP2 method for the Ace-(Ala) n -Nme system in vacuo. [14] It is not straightforward to extract the H-bond energy as a part of a large molecule, where the donor and acceptor atoms are linked through several covalent bonds. In fact, 12 covalent bonds link the acceptor atom, O, and the donor atom, H, in the backbone α-helical H-bond between the i-th and (i+4)-th amino-acid residues. Namely, the backbone α-helical H-bond is not a simple, additive pair-wise interaction, because it includes many-body effects with nonadditive natures.
For that purpose, we employed the MTA developed by Deshmukh and Gadre, who showed that it is possible to estimate the intramolecular backbone H-bond energies of 3 10 -helices in several model polypeptides. [16] Here, we used MTA to systematically compute the backbone H-bond energies in several α-helical peptide models.
The total energy of Ace-(Ala) n -Nme is approximated by dividing-and-conquering the energies of the fragments by MTA, using the following equation: In eq. (1), each sum was taken for the possible combinations for (a, b, …), where 1 ≤ a, b, … ≤ (n + 1). An example (n = 3) of the fragment models is shown in the Supporting Information Figure S1. The total energies of the systems are well approximated by the combinations of all possible fragments, as described below, and thus we computed the energy of the entire system (G 0 in Fig. 2), instead of the combination in the original MTA method. [16] We also computed the energy of a peptide fragment lacking the acceptor group (G 1 ), the energy of a fragment lacking the donor group (G 2 ), and the energy of a fragment lacking both the acceptor and donor groups (G 12 ), as shown in Figure 2. The H-bond energy, E HB , and the electron density change upon H-bond formation, Δρ MTA , were estimated as follows: We also calculated the stabilization energy (SE) in eq. (4), which we introduced in our previous paper, [14] by using the same DFT B97D/6-31+G(d) method: where E MH N&C is the total energy of the MH model with both the N-and C-terminal N-methylacetamide molecules (Fig. 1B). E MH N and E MH C are the total energies of the two separated N-methylacetamide molecules at the N-and C-termini, respectively. Namely, ΔE SE total is the energy difference between the total energy of the MH model with the α-helical H-bond and that where the H-bond acceptor and the donor are separated infinitely. All the SE values were corrected for the basis set superposition error (BSSE) by the counterpoise method of Boys and Bernardi. [18] The ordinary electron density change upon H-bond formation, Δρ SE , was computed as follows: As references, we also computed the MM interaction energies for the corresponding H-bond interactions, using the following equation (6), where i and j are four atoms attributed to peptides I and J, respectively: {C, O, N, and H}. A ij and B ij are the Lennard-Jones coefficients, r ij is the distance between i-th and j-th atoms, and q i is the atomic partial charge of the i-th atom. Here, AMBER ff99SB force-field parameters [15] were used. Electron density changes were computed using the cube files in the Gaussian09 program package, [17] and the figures of the molecules with the electron density changes were produced by UCSF Chimera. [19] Results

H-bond energies in α-helices
The total energies of Ace-(Ala) n -Nme estimated by MTA (E MTA ), computed by eq. (1), coincided well with the ordinary total energies E(F 0 ) of the complete AH models. In fact, the differences in the values calculated by E(F 0 ) and MTA, E MTA -E(F 0 ), are indicated in Table 1, and all of them were <0.09 kcal/mol, representing about 0.00001% of the total energies. These differences are similar to that obtained in the previous study of the 3 10 -helix, where the difference was 0.11 kcal/mol. [16] In Figure 3A, the H-bond energies between the backbone donors and acceptors are plotted for individual pairs for the models AH3 to AH8, depending on the distance between the donors and the acceptors for the AH models, ST models, and MH models by the MTA method, and the classical H-bond energies given by the MM calculation. The colors indicate the α-helical structures (Ace-(Ala) n -Nme) with different lengths (3 ≤ n ≤ 8), as indicated in the caption of Figure 3. The structural variations were caused by the energy minimization procedures for the entire α-helical conformations, as mentioned in the Methods section. The actual energy values are summarized in Table S1 in the Supporting Information. During the energy minimization procedure, the distance between O atom of C O at i-th residue and H atom of HN at (i+4)-th residue became longer than 1.98 Å in the initial structure, and this structure modification was slightly shorter in the longer α-helices.
As shown in Figure 3A and in the Supporting Information Table S1, the N-terminal helical turns were largely deformed from the initial conformations providing different H-bond lengths during the energy minimization procedure. Moreover, the angle between the two vectors of the carbon atom to the oxygen of the carbonyl group of the i-th residue and the nitrogen to the hydrogen of the amide group of the (i+4)-th residue was also deformed from the initial value at the N-terminus, as shown by the open and thin symbols in Figure 3B. Those vector pairs formed angles larger than 170 , and most of them were close to 180 . At the third H-bond, this angle was slightly smaller, by about 5 , than those of many other typical vector pairs, which were about 168 (Fig. 3B).
Here, the H-bond energies in all of the models correlated well. In fact, when the H-bond energies calculated by the MTA method for the AH and ST models are plotted against those for the MH model, the H-bond energies calculated for the AH and ST models strongly correlated with that for the MH model, as shown in Figure 4, with Pearson's correlation coefficients of 0.885 and 0.987, respectively.
The MM values have also a high correlation coefficient, 0.975, with the MH model, and even the absolute MM values are very close to the H-bond energies calculated by the MH models. In contrast, the H-bond energies obtained by both the AH and ST models remarkably deviated from those calculated by the MH models.
The reaction field due to high-dielectric water generally tends to shield the electrostatic interactions. In order to examine the actual shielding effect to the H-bonds in α-helix by water, the polarizing continuum method (PCM) with the high dielectric constant (ε = 78.3553) was simply applied to the MTA method using the same DFT function by Gaussian09. [20] As shown in the Supporting Information Table S4 and Figure S6, the H-bond energies became about 1 kcal/mol weaker than those in vacuo for the current α-helical structures.

SEs in MH models
The SE values defined in eq. (4) for the MH models, ΔE SE total , are also plotted in Figure 4, and the actual values are provided in Table S1. The SE values correlated very well with the corresponding H-bond energies in the MH models, with a Pearson's correlation coefficient of 0.995. The SE values were 0.93 kcal/mol lower than the corresponding H-bond energies in the MH models except for the N-termini, where the N-terminal backbone structures were largely deformed and their SE values were 0.86 kcal/mol lower than the H-bond energies in the MH models.

Electronic structures around the H-bond donors and acceptors
In addition to the H-bond energies, the MTA method can approximate the electronic structures around the H-bond Table 1. Total energies by the original DFT computation, E(F 0 ), and the MTA method (E MTA ) for the AH3, 4, 5, 6, 7, and 8 structures, corresponding to Ace-(Ala) n -Nme with n from 3 to 8. The differences are also shown.  alanine octamer, Ace-(Ala) 8 -Nme, was provided by eq. (5), and it is shown in Figure 5A. The corresponding electron density changes for the AH, MH, and ST models, Δρ MTA , computed by eq. (3), are also shown in Figures 5B-5D, respectively. Here, yellow and magenta colors show the negative and positive contour surfaces, respectively. It is obvious that the Δρ MTA values are all similar to the Δρ SE values, where the electron density increases around the oxygen atom of the carbonyl (C O) group at the i-th residue and decreases around the hydrogen atom of the amide (N─H) group at the (i+4)-th residue. For the other structures from  to , the Δρ MTA values of the AH models are shown in the Supporting Information Figures S2A-S2E. The differences in the electron density changes between the AH and MH models and those between the ST and MH models were further computed, respectively: ΔΔρ In Figures 5E and 5F, the ΔΔρ

Discussion H-bond energies by the MTA method
The energies of Ace-(Ala) n -Nme provided by the MTA method are approximated values. However, as shown in Table 1, the differences in the values between the ordinary total energies for the complete F 0 and those obtained by MTA, E MTA -E(F 0 ), are very small. They are also small even in comparison to the Hbond energies, which are about −3 to -4 kcal/mol in this study. Thus, the H-bond energies estimated by the current MTA method should be quantitatively reliable, with about 3% errors, for discussing the H-bond interactions in α-helices.
The accuracy of the MTA method largely relies on the borders separating molecular segments. Historically, this issue was recognized as the "nearsightedness of electronic matter (NEM)" [21] to divide-and-conquer large molecular systems in general. Using the theoretical computations on the basis of the linear response function, the sp 3 junction was the most suitable location for partitioning peptide systems. [22,23] In the current MTA procedures, all the fragmentations followed the sp 3 junction mechanism to block the propagation of the electron density deviation.
In the MH models where the two peptide groups for hydrogen donors and acceptors are separated without any covalent bonds, the H-bond energies should directly correspond to the SEs including the BSSE corrections. [18] In fact, an almost perfect correlation appeared between the H-bond energies and the SEs, as shown in Figure 4 and in the Supporting Information Table S1, and the differences were always the same, 0.93 kcal/mol, except for the N-termini where the backbone structures were largely deformed during the energy optimization procedures. Those differences are due to other interactions among the methyl groups that capped the N-and C-terminal peptide groups than the H-bond energies given by the current MTA analysis using the MH models. In addition, as shown in Figures 5A and 5C, the change in the electron density upon Hbond formation, Δρ SE in eq. (5), is also well approximated by Δρ MTA MH in eq. (3).
Thus, the current MTA method provides good approximations of the H-bond energies and electronic structures in MH models, and so it is expected to give a reliable analysis for the H-bond interactions of the AH and ST models as well, as shown in Figures 5B and 5D. Morozov et al. [24] reported a similar approach to analyze the cooperativity of α-helix formation by using separated α-helical peptide fragments. In particular, their model including the short-range contribution was, in principle, designed to compute the dimerization energies using the SEs. Namely, their shortrange interactions did not correctly account for the nonadditive many-body interactions and ignored the effects of the α-helical backbone atoms linking the H-bond acceptor and donor.

Analysis of the electronic structures
It is clear from Figure 4 and from the Supporting Information Table S1 that the H-bond energies for the AH models have similar values to those for the ST models, although the former ones tend to be slightly weaker than the latter ones. In contrast, the H-bond energies in the AH and ST models significantly deviated from those in the MH models, although the Pearson's correlation coefficients were very high. As shown in Figures 5E and 5Fand in the Supporting Information Figures S3 and S4, the electronic structures around the H-bonds in the AH and ST models distinctively deviated from those in the MH models, with the depolarization of the hydrogen donor and acceptor groups.
Thus, the phenomenon should be caused by the precise electronic structures in the ST and AH models. In order to analyze the origin of this phenomenon, we focused on the six model structures of AH8-1 to AH8-6. We found that the distances between the oxygen atoms in the carbonyl group of the i-th and (i+1)-th residues in the six H-bond pairs are short, 3.510AE0.144 Å. In addition, those between the hydrogen atoms in the amide group of the (i+3)-th and (i+4)-th residues are also short, 2.676AE0.038 Å. These short distances suggest that the carbonyl oxygen of the i-th residue has less electron density, and the amide hydrogen of (i+4)-th residue has more electron density, as revealed in Figures 5E and 5F.
By applying the Hirshfeld population analysis, [25][26][27] the electronic structures in the six ST models were analyzed around the carbonyl oxygen of the i-th residue and the amide hydrogen of the (i+4)-th residue, in comparison to those in the MH models, in which no neighboring carbonyl (C O) or amide (NH) groups exist. As shown in the Supporting Information Table S2C, for the G 2 fragment of the ST model (Fig. 2C), which lacked the H-bond donor group at the (i+4)-th residue but included the effect of the C O group of the (i+1)-th residue, the Hirshfeld atomic charge of the carbonyl oxygen of the i-th residue was 0.0303eAE0.0036e larger than that in the MH model. The contribution of the amide hydrogen of the (i+4)-th residue for the G 1 fragment of the ST model, which lacked the H-bond acceptor group at the i-th residue but included the effect of the amide group of the (i+3)-th residue, was 0.0154eAE0.0010e less than that in the MH model. Similar Hirshfeld charge changes were also observed for the entire G 0 systems, where H-bonds are formed between the C O group of the i-th residue and the NH group of the (i+4)-th residue, as shown in the Supporting Information Table S2A.
These depolarization effects also correlate with the local dipole moments. The local inter-atomic dipole moment for a system composed of N charges {q k | k = 1, …, N} is described by eq. (9) where <q > is the average of the Hirshfeld atomic charges and < r ! > is their center position. For the C O group of the i-th residue and the NH group of the (i+4)-th residue, the local dipole moments become simple, as shown in eqs. (10) and (11): Here, q k C and q k O are the Hirshfeld atomic charges of the C and O atoms in the C O group of the k-th residue, and q k H and q k N are those of the NH group of the k-th residue, respectively. HN |, are used for the following discussion. Due to the closely located electric dipole μ i+1 CO at the (i+1)-th residue to μ i CO at the i-th residue in the parallel direction, when an α-helical conformation is formed, the dipole moments should decrease by their repulsive interaction. In the same way, due to the closely located electric dipole μ i+3 HN at the (i+3)-th residue to μ i+4 HN at the (i+4)-th residue in the parallel direction in an α-helix, those dipole moments should also decrease. From the Hirshfeld population analysis, as shown in Table 2B, the average ratio of μ i CO in the ST model to that in the MH model was 0.938AE0.002, and the average ratio of μ i+4 HN in the ST model to that in the MH model was 0.917AE0.002, when the H-bonds were not formed between the C O group of the i-th residue and the NH group of the (i+4)-th residue. As shown in Table 2A, when Hbonds were formed between the C O group of the i-th residue and the NH group of the (i+4)-th residue, the average ratio of μ i CO in the ST model to that in the MH model was 0.944AE0.002, and the average ratio of μ i+4 HN in the ST model to that in the MH model was 0.946AE0.003. Namely, these dipole moments were reduced by about 5-6 and 5-8%, respectively, depending on the backbone α-helical conformation. Although such population analyses may include some ambiguities for the absolute values of the atomic charges and dipole moments, these tendencies accurately reflect the changes of the electronic structures shown in Figures 5F and in the Supporting Information Figure S4, as the depolarization effect.
In order to examine the above phenomena, two other models, the HT N (N-terminal half-turn) and HT C (C-terminal halfturn) models, were constructed. In the HT N model, the (i+2)-th and (i+3)-th Ala residues were both deleted from the ST model and capped by methyl groups, as shown in the Supporting Information Figure S5A. In the HT C model, the (i+1)-th and (i+2)-th Ala residues were both deleted, as shown in the Supporting Information Figure S5B. Namely, the carbonyl group of the (i+1)-th residue is included in the HT N model, and the amide group of the (i+3)-th residue is included in the HT C model. Consequently, the H-bond energies for the HT N and HT C models were located in the middle between the corresponding ST and MH models, as shown in Figure 6 by the filled and open triangles, respectively. The actual energy values are summarized in the Supporting Information Table S3. A careful investigation of the H-bond energies of the AH and ST models in the Supporting Information Table S1 and Figure 6 revealed that the H-bond energies in the AH models are always slightly weaker than those in the ST models, except for the Ntermini or C-termini. These effects can be caused by the successive carbonyl group of the (i-1)-th residue and the amide group of the (i+5)-th residue, in a similar manner to the effects of the successive carbonyl group of the (i+1)-th residue and the amide group of the (i+3)-th residue in the opposite directions. In fact, the distances between the oxygen atoms in the C O groups of the (i-1)-th and i-th residues in the five H-bond pairs are 3.506AE0.037 (Å), and those between the hydrogen atoms in the NH groups of the (i+4)-th and (i+5)-th residues are 2.714AE0.038 (Å), as shown in the Supporting Information Table S2B.
Two additional models, the AP N (N-terminal additional peptide) and AP C (C-terminal additional peptide) models, were also considered, where an Ace-Ala group and an Ala-Nme group were added to the N-and C-terminus of the MH model, respectively (Supporting Information Figs. S5C and S5D). The AP N model has an interaction between the successive C O groups at the (i-1)-th and i-th residues, and the AP C model has another interaction between the successive NH groups at the (i+4)-th and (i+5)-th residues. Their H-bond energies computed by the MTA method are shown in Figure 6 and in the Supporting Information Table S3. Both of them have the middle H-bond energies between the corresponding ST and MH models, suggesting another depolarization effect.
These phenomena were also analyzed with the local electric dipole μ i CO at the i-th residue in the AP N model, and μ i+4 HN at the (i+4)-th residue in the AP C model. Consequently, the ratio of μ i CO in the AP N model to that in the MH model was 0.971AE0.004, and the ratio of μ i+4 HN in the AP C model to that in the MH model was 0.977AE0.001, when the H-bonds were not formed between the C O group of the i-th residue and the NH group of the (i+4)-th residue (Table 2B). When H-bonds were formed between the C O group of the i-th residue and the NH group of the (i+4)-th residue, the ratio of μ i CO in the AP N model to that in the MH model was 0.973AE0.004, and the ratio of μ i+4 HN in the AP C model to that in the MH model was 0.985AE0.001 (Table 2A). These dipole moments were both also reduced by about 2-3%, and their contributions to the total Hbond energies are smaller than those of μ i+1 CO or μ i+3 HN . As shown in Figure 6 and in the Supporting Information Table S3, the H-bond energies due to the surrounding carbonyl groups around the i-th residue and those of the amide group around the (i + 4)-th residue were not additive. Namely, the simple summations of the energy differences of the H-bond energies between the HT N and MH models and between the HT C and MH models are always larger than the direct differences between the ST and MH models. Similarly, simple summations of the differences in the H-bond energies between the AP N and MH models and between the AP C and MH models could also overestimate the differences between the AH and ST models. Thus, they should be considered as nonadditive many-body effects.
Finally, the putative effects of the helical dipoles at the backbone peptide planes that are far from the target H-bond locations were also investigated. Although the neighboring peptide dipoles at the (i-1)-th and (i+5)-th residues slightly contributed to the H-bond energy of the unit α-helix from the i-th to (i+4)-th residues, as discussed earlier, there was no significant dependence of the H-bond energies on the helix length, as shown in Figure 4 and in the Supporting Information Table S1. The farther helical dipoles do not seem to affect the electronic structures of the H-bond donors and acceptors, although they contribute to the cooperative nature of the α-helix formation. [24] Thus, we can conclude that the H-bond energies of the α-helix, as in the AH and ST models, are generally weaker than those of the separated H-bonds, as in the MH model, due to the depolarized electronic structures around the carbonyl oxygen of the i-th residue and the amide hydrogen of the (i+4)-th residue. Such depolarizations redistribute the electron density and are caused by the local electronic interactions in their neighborhood inside the α-helical structure. Similar H-bond energy changes depending on the peptide backbone structures were also found in the antiparallel β-sheet models in our previous paper [14] and others. [28,29] When the SEs were computed, the odd-numbered β-sheet models had weaker SEs by forming smaller hydrogen bond ring structures, and the even-numbered β-sheet models had stronger SEs by forming larger hydrogen bond ring structures. [14] Toward improvement of the H-bond energy by the classical force field The current analysis of the depolarization of the electronic structure at the carbonyl group of the i-th residue and the amide group of the (i+4)-th residue in an α-helix revealed the