Quantum chemical studies on hydrogen bonds in helical secondary structures

We present a brief review of our recent computational studies of hydrogen bonds (H-bonds) in helical secondary structures of proteins, α-helix and 310-helix, using a Negative Fragmentation Approach with density functional theory. We found that the depolarized electronic structures of the carbonyl oxygen of the ith residue and the amide hydrogen of the (i + 4)th residue cause weaker H-bond in an α-helix than in an isolated H-bond. Our calculations showed that the H-bond energies in the 310-helix were also weaker than those of the isolated H-bonds. In the 310-helices, the adjacent N–H group at the (i + 1)th residue was closer to the C=O group of the H-bond pair than the adjacent C=O group in the 310-helices, whereas the adjacent C=O group at the (i + 1)th residue was close to the H-bond acceptor in α-helices. Therefore, the destabilization of the H-bond is attributed to the depolarization caused by the adjacent residue of the helical backbone connecting the H-bond donor and acceptor. The differences in the change in electron density revealed that such depolarizations were caused by the local electronic interactions in their neighborhood inside the helical structure and redistributed the electron density. We also present the improvements in the force field of classical molecular simulation, based on our findings. Supplementary Information The online version contains supplementary material available at 10.1007/s12551-022-01034-5.


Introduction
Proteins are found in all living organisms and are involved in almost all biological activities such as catalysis, molecular recognition, and material transport (Liljas et al. 2017;Branden and Tooze 2012). Since protein functions are strongly correlated with their three-dimensional structures, understanding the three-dimensional structure of proteins and their dynamic behaviors is essential for various scientific fields including chemistry, biology, medicine, agriculture, and food industry.
Proteins are biopolymers consisting of a large number of amino acids held together by peptide bonds. Protein structures are hierarchical, with distinct levels of structures (Holde Van et al. 1998), which represent increasing levels of complexity and include primary, secondary, tertiary, and quaternary structures. Secondary structure is the local and regular structure of a protein, including α-and 3 10 -helices, β-strands, and β-and γ-turns. These secondary structures are properly assembled to form tertiary structures. Thus, secondary structures can be considered the building elements of protein structures.
In secondary structures, a helix (e.g., α-helix and 3 10 -helix) is a group of residues that repeatedly rotate and rise along an axis. It is the most observed secondary structure. We can classify it into different helical conformations. The α-helix is found in 31% of all secondary structures and is the most widely recognized helical structural element in fibrous and globular proteins (Barlow and Thornton 1988).
The second most common helical structure is the 3 10 -helix, which occupies 4% of the secondary structures.
These helices are stabilized by hydrogen bonds (H-bonds) formed between the amide hydrogens (the H-bond donors) and the carbonyl oxygens (the H-bond acceptors) of peptide bonds. H-bonds are one of the most important noncovalent interactions for chemical and biological phenomena (Saleh et al. 2012;Tantardini 2019;Tantardini et al. 2020). The significance of H-bonds in the secondary structures of proteins was recognized early (Pauling et al. 1951;Eisenberg 2003). The α-helix has 3.6 residues per turn and is a right-handed helix. In the α-helix, a H-bond is formed between the peptide carbonyl group at residue i and the peptide amino group at residue i + 4. In contrast, the 3 10 -helix has three residues per turn and is a right-handed helix. It has a H-bond between the peptide carbonyl group at residue i and the peptide amino group at residue i + 3, resulting in a tighter packing of the backbone compared with the α-helix (Fig. S1). Many computational chemists have studied H-bonds in secondary structures at various levels of theoretical depth. Wieczorek and Dannenberg (2003a, b) investigated H-bond cooperativity and the energetics of α-helices, suggesting that various factors contribute to their stability. Morozov et al. (2006) evaluated the origin of cooperativity in forming α-helices. Wu and Zhao (2001) studied the role of cooperativity in the formation of α-helices by performing theoretical calculations on α-helix models constructed using a simple repeating unit method. Parthasarathi et al. (2007) studied H-bond interactions in an α-helix model using the atom-in-molecules method. Ismer et al. (2008) investigated the temperature dependence of the stability of α-, π-, and 3 10 -helices compared with a fully extended structure using density functional theory (DFT) and harmonic approximation. However, a simple physicochemical theory accounting for helical secondary structural features of proteins is still immature.
An accurate and quantitative evaluation of H-bonds is also important for molecular dynamics (MD) simulations to investigate the dynamical behavior and folding process of proteins. Historically, it has been noticed since many years ago that each force field used in classical MD simulations shows a specific tendency to form an α-helix or a β-strand. For example, the AMBER C96 force field, which was developed just after the original AMBER force field , preferred extended structures contrast to the α-helical preference of the latter one as mentioned by Kollman et al. (1997). This phenomenon has been repeatedly reported by many authors (Sakae and Okamoto 2003;Yoda et al. 2004a, b;Best et al. 2008;Best and Hummer 2009;Piana et al. 2011). Usually, this preference for force fields on the secondary structure formation is not a significant problem in the MD simulations of rigid globular protein structures. However, it has become a critical issue in understanding functionally important conformational changes (Higo et al. 2011;Shirai et al. 2014;Chebaro et al. 2015;Nishigami et al. 2016) in the folding simulations of flexible disordered regions (Higo et al. 2011;Chebaro et al. 2015) and long loops between secondary structures (Shirai et al. 2014;Nishigami et al. 2016). Yoda et al. (2004a, b) performed MD simulations to compare the secondary structural properties of commonly used force fields, finding that MD simulations with the AMBER ff94  and ff99 (Wang et al. 2000) force fields were in remarkable agreement with experimental data for α-helical polypeptides but not for β-hairpin polypeptides. Numerous attempts have been made to overcome this problem, such as increasing the torsional energies, rearrangements (Kamiya et al. 2005;Buck et al. 2006;Fujitani et al. 2009;Robustelli et al. 2018), and developing polarized charge models (Patel and Brooks 2004;Lopes et al. 2009). Regardless, the reasons behind the use of these methods remain unclear, and elucidation requires understanding the energy of hydrogen bonding in the secondary structure.
Here, we briefly review our computational studies of H-bonds in helical secondary structures of proteins, α-helix and 3 10 -helix, using a Negative Fragmentation Approach (NFA) with DFT (Kondo et al. 2019(Kondo et al. , 2022. We will also discuss the modification of the force field, approximating the H-bond energies, revealed by our findings.

Model constructions
We constructed whole-helix (WH) models of an α-helix and a 3 10 -helix, denoted as WH alpha -n and WH 3_10 -n models, respectively. In the model construction, we used poly-alanine amino acids capped with an acetyl group (ACE) and an N-methyl amide group (NME), denoted as ACE-(Ala) n -NME (n = 2-7 for the 3 10 -helices and n = 3-8 for the α-helices). The backbone torsion angles, φ and ψ, of the WH alpha -n and WH 3_10 -n models were set to their ideal values, as described in biochemistry textbooks, namely, φ = − 57° and ψ = − 47° for the WH alpha -n models and φ = − 49° and ψ = − 26° for the WH 3_10 -n models (Arnott and Dover 1967;Petsko and Ringe 2004;Kuster et al. 2015). These structures were optimized in the gas phase by fixing the backbone dihedral angles at the aforementioned values and minimizing the total electronic energies. One to six backbone H-bonds exist in the optimized WH alpha -n and WH 3_10 -n models. The sth H-bond in these models, counting from the N-terminus, is represented by n-s.
To understand the characteristics of the H-bond energies in α-helices and 3 10 -helices, we constructed two types of simplified models. One is a single-turn (ST) model, denoted as the ST alpha -n and ST 3_10 -n model composed of ACE-(Ala) 3 -NME and ACE-(Ala) 2 -NME, respectively. In this model, the H-bond donor and acceptor are linked with the helical backbone atoms. The other is a minimal H-bond (MH) model, denoted as MH alpha -n and MH 3_10 -n models comprising two separated N-methyl acetamide molecules and mimicking a single H-bond between the C=O and N-H groups in the backbone. In the MH alpha -n and MH 3_10 -n models, the two peptide groups forming a H-bond, hydrogen donors and acceptors, were separated without linking the helical backbone atoms. The atomic positions of these simplified models were the same as those of the corresponding WH alpha and WH 3_10 models, except for the N-and C-terminal capping groups. Figure 1a shows the molecular structures of the WH alpha -5, ST alpha -5, and MH alpha -5 models for an α-helix and those of the WH 3_10 -4, ST 3_10 -4, and MH 3_10 -4 models for a 3 10 -helix as examples. The individual H-bond energy for these models was calculated in the same manner as that for each backbone H-bond in the WH models, as described below. The H-bond energies of these models were then compared to each other.

Validity of DFT exchange-correlation functionals for the calculation of H-bond energy of secondary structure
In recent years, DFT has become accepted as an alternative approach for the post Hartree-Fock (HF) methods such as Møller-Plesset perturbation theory (Head-Gordon et al. 1988) and coupled cluster theory (Scuseria and Schaefer 1989). In previous studies (Takano et al. 2011(Takano et al. , 2012, we showed the importance of assessing the validity of Fig. 1 a Structures of the WH alpha -5, ST alpha -5, and MH alpha -5 models for an α-helix and those of the WH 3_10 -4, ST 3_10 -4, and MH 3_10 -4 models for a 3 10 -helix as examples. b Fragment structures and schematic picture for calculations of H-bond energy of 3-1 of the WH alpha -3 model with NFA various DFT exchange-correlation functionals. The DFT exchange-correlation functionals was also validated for the H-bond energies of the ACE-(Ala) n -NME system. We chose the B97D (B97 functional with Grimme's D2 dispersion schemes) exchange-correlation functional (Grimme 2006) with 6-31+G(d) basis sets because it provided the H-bond energies of an ACE-(Ala) n -NME dimer comparable with the MP2 method in the calculation of the H-bond interaction energies of the ACE-(Ala) n -NME system in the gas phase (Takano et al. 2016). This method was applied to H-bond energy calculations of the helical secondary structures composed of ACE-(Ala) n -NME.

Negative fragmentation approach
Since it is not straightforward to calculate a H-bond energy of a molecule, where the donor and acceptor atoms are linked through covalent bonds, we utilized the NFA, which is the modified version of the Molecular Tailoring Approach (MTA) developed by Deshmukh and Gadre (2009). As shown in Fig. 1b, each H-bond energy, E HB , in ACE-(Ala) n -NME is calculated by the following equation in the NFA: E sys , E noA , E noD , and E noAD represent the total electronic energy of the entire system, the system lacking the acceptor group, the system without the donor group, and the system lacking both the acceptor and donor groups, respectively.
In the NFA, we used the total electronic energy of the entire system E Total as E sys (Kondo et al. 2019(Kondo et al. , 2022, while the energy of the entire system was estimated using the energies of all fragments in the original MTA (Deshmukh and Gadre 2009). The total energies of the WH alpha -n models estimated by MTA E MTA coincided well with the E Total values of these models. The differences in the calculated values of E Total and E MTA , E MTA − E Total , were less than 0.09 kcal/ mol. These differences are similar to that obtained in the previous study (Deshmukh and Gadre 2009) for the 3 10 -helix (0.11 kcal/mol).
In addition to the H-bond energies, the NFA can approximately represent the change of electronic structures upon H-bond formation. The change in electron density upon H-bond formation, Δ HB , was evaluated as follows (Kondo et al. 2019(Kondo et al. , 2022: To examine the difference in the H-bond energy between the WH and MH models and between the ST and MH models in the context of their electronic structures, the (2) Δ HB = sys − noA − noD + noAD differences in the change in electron density were computed using Eqs. 3 and 4, respectively: As an advantage of NFA, no modifications of program code are required, though multiple calculations are needed. It indicates that we can utilize the NFA with any quantum chemical calculation programs.
For comparison and improvement of the force field, we also computed the H-bond energies based on the molecular mechanics (MM) with the AMBER ff99SB force field parameters (Wang et al. 2000), E MM HB , for the corresponding H-bonds as in the following equation: where i and j are the atoms constituting the peptide group of an acceptor and a donor of an H-bond, respectively: {C, O, N, and H}. A ij and B ij are the Lennard-Jones coefficients, r ij is the distance between atoms i and j, and q i is the atomic partial charge of the atom i.

H-bond energies in helical model systems
In order to understand the effects of the helical secondary structures on the H-bond energies, we compared the H-bond energies for the WH alpha -n, ST alpha -n, and MH alpha -n models of α-helices and the WH 3_10 -n, ST 3_10 -n, and MH 3_10 -n models of the 3 10 -helices. Since the H-bond energies strongly depend on the spatial arrangement of the H-bond donor and acceptor atoms, we plotted the H-bond energies of the WH and ST models against those of the MH models for the αand 3 10 -helices in Fig. 2a and b, respectively, to cancel out the effect of the orientations of H-bond donor and acceptor (Kondo et al. 2022).
The H-bond energies obtained by the WH-n and ST-n models, E WH HB and E ST HB , remarkably deviated from those calculated by the MH models. In the α-helices, the ST alpha -n model reproduced the H-bond energies of the WH alpha -n model. In contrast, the MH alpha -n model provided more stable H-bond energy than the WH alpha -n model (Fig. 2a), indicating that the adjacent residues covalently connecting the H-bond donor and acceptor destabilized the H-bond in the WH alpha -n model. In the 3 10 -helices, the ST 3_10 -n models also destabilized the H-bond as well as the α-helices compared to the MH 3_10 -n models but failed to provide the ij equivalent H-bond energies of the WH 3_10 -n. In particular, the H-bond pairs adjacent to the N-or C-terminal H-bond of the WH 3_10 -n models were strongly destabilized compared to those of the ST 3_10 -n models, resulting in the worse correlation of the WH 3_10 -n and ST 3_10 -n models. Details are discussed in our previous study (Kondo et al. 2022). It suggests , for 5-2 of the WH alpha -5 model and for 4-2 of the WH 3_10 -4 model. The yellow surfaces represent the contour surfaces at −0.001 au, and the magenta ones are those at +0.001 au. The atoms in the whole WH models are shown by green wire and those in the MH models are shown using the stick model with CPK colors. The difference in the change in electron density between the WH and ST models, ΔΔ  . Figures d, f, and h were slightly modified from the original figures, which appeared in the previous paper by Kondo et al. (2022) that the destabilization of the H-bond of the 3 10 -helices is partly due to the helical backbone atoms linking the H-bond donor and acceptor. However, there should also be other factors leading to unstable H-bond energies.
We compared QM ( E WH HB ) and MM ( E MM HB ) calculations for the H-bond energies of the WH alpha -n and WH 3_10 -n models (Fig. 2a, b). In the WH alpha -n models, the E WH HB values strongly correlated with the E MM HB values, but the MM calculations overestimated the stability of the H-bond energies by ~ 1 kcal/mol. The E MM HB values of the WH alpha -n models almost coincided with the H-bond energies calculated by the MH models, E MH HB . This is because the current force field parameters of amino acids are adjusted based on the amino acid monomer. In contrast, the H-bond energies evaluated with the QM calculations were closer to those with the MM calculations for the WH 3_10 -n models than for the WH alpha -n models. However, the correlation between the QM and MM calculations was weak, unlike the WH alpha -n models. In contrast to α-helices, the E MH HB values were more stable than the E MH HB values in 3 10 -helices. This is because the H-bonds in the MH 3_10 -n are much shorter than those in the MH alpha -n models, thus having stronger quantum nature that the classical force field cannot describe.

Electronic structures around the H-bond donors and acceptors
The electron density changes for the WH alpha -5 and WH 3_10 -5 models, Δ WH alpha HB and Δ WH 3_10 HB , were calculated with Eq. 2 and are shown in Fig. 2c and d, respectively. Here, yellow and magenta colors show the negative and positive contour surfaces, respectively. From the Δ WH alpha HB values, the electron density increased around the oxygen atom of the C=O group at the ith residue and decreased around the hydrogen atom of the N-H group at the (i + 4)th residue, as shown in Fig. 2c. The Δ WH 3_10 HB values showed that the electron density increased in the vicinity of the oxygen atom of the C=O group at the ith residue and that it decreased in the vicinity of the hydrogen atom of the N-H group at the (i + 3)th residue, as illustrated in Fig. 2d In Fig. 2e and g, the differences in electron density change between the WH alpha -5 and MH alpha -5 models and between the ST alpha -5 and MH alpha -5 models, ΔΔ WH alpha −MH alpha HB and ΔΔ ST alpha −MH alpha HB , respectively, were shown for the α-helical turn 5-2. Here, green and orange colors show the negative and positive contour surfaces, respectively. The electron density changes near the oxygen atom of the C=O group at the ith residue in both the WH alpha and ST alpha models were smaller than that in the MH alpha model, implying that the negative polarization of the oxygen atom was weakened by the backbone atoms linking the H-bond donor and acceptor. In contrast, the electron density near the hydrogen atom of the N-H group at the (i + 4)th residue increased in the WH alpha and ST alpha models, as compared with that in the MH alpha model. In addition, the ΔΔ WH alpha −MH alpha HB value was similar to the ΔΔ ST alpha −MH alpha HB values. It indicates that the weaker positive polarization of the hydrogen atom is mainly due to the adjacent residues connecting the H-bond donor and acceptor. We found that the distances between the oxygen atoms in the carbonyl group of the ith and (i + 1)th residues in the H-bond pairs were short (3.510 ± 0.144 Å). In addition, those between the hydrogen atoms in the amide group of the (i + 3)th and (i + 4)th residues were also short (2.676 ± 0.038 Å). These short distances caused the depolarization of both the carbonyl oxygen of the ith residue and the amide hydrogen of the (i + 4)th residue, as revealed in Fig. 2e and g. The depolarized electronic structures around the carbonyl oxygen of the ith residue and the amide hydrogen of the (i + 4)th residue generally resulted in weaker H-bond energies for the α-helix, as in the WH alpha and ST alpha models, than for the separated H-bonds, as in the MH alpha model. Such depolarizations redistributing the electron density were caused by the local electronic interactions in their neighborhood inside the α-helical structure. values, implying other factors besides the helical backbone atoms linking the H-bond pairs that caused the depolarization of the H-bond donor and acceptor in the 4-2 pair. We now discuss why the H-bonds in the ST 3_10 -n model were destabilized in comparison to those in the MH 3_10 -n model. We found that the C=O group at the ith residue and the N-H group at the (i + 3)th residue of the H-bond were depolarized in the ST 3_10 -n model, as shown in Fig. 2h. This depolarization could be caused by the helical backbone atoms linking the H-bond pair. Whereas the C=O group at the (i + 1)th residue was involved in depolarization in α-helices, the C=O group of the H-bond pair was closer to the adjacent N-H group at the (i + 2)th residue (2.799 ± 0.033 Å), than to the adjacent C=O group at the (i + 1)th residue (3.442 ± 0.025 Å) in the 3 10 -helices. Therefore, in the 3 10 -helices, the adjacent N-H group may cause the depolarization of the H-bond acceptor, resulting in only little destabilization of the H-bond.

Toward improvement of the H-bond energy by the classical force field
Our calculations for α-helices and 3 10 -helices with the NFA revealed that their H-bond energies are affected by depolarization and polarization due to the local dipole of the neighboring backbone peptide groups. Based on our results, we constructed a model for the classical force field, in which the atomic partial charges ( q i N , q i H , q i C , and q i O ) of the N-H and C=O groups of the ith peptide group were not constant but were changed by the neighboring peptide groups, respectively.
where q 0 N , q 0 H , q 0 C , and q 0 O are the original atomic partial charge, and the amounts of changes in the atomic partial charges of the ith peptide group; i N , i H , i C , and i O , are functions depending on the backbone structure of the neighboring peptide group. This change may be represented by the interaction energy of the classical mechanics, U � ⃗ i X , � ⃗ j Y between the ith backbone dipole, � ⃗ i X (X = NH or CO), and the neighboring backbone dipole, � ⃗ j Y (Y = NH or CO) ( j ≠ i ). To avoid double counting interaction energies already taken in contributions in the restrained electrostatic potential (RESP) approach (Bayly et al. 1993;Cieplak et al. 1995), the following formula was used to obtain the change in the original atomic partial charge.
where HO , HH , OO , and OH are fitting parameters. To ensure that the overall charge does not change, local neutrality conditions were imposed, and i N and i C are defined as follows.
We determined the parameters HO , HH , OO , and OH , by using the H-bond energies calculated with NFA, HB , for both α-helices and 3 10 -helices, and those of the AMBER ff99SB force field parameters (Wang et al. 2000 was computed from the modified charges, q N , q H , q C , and q O , given by Eqs. (6)-(9) as the changes from the original charges, q 0 N , q 0 H , q 0 C , and q 0 O , which were −0.4157, 0.2719, 0.5973, and −0.5679, respectively, taken from the AMBER ff99SB force field (Wang et al. 2000). The parameters modifying charges, i N , i H , i C , and i O , in Eqs. (6)-(9) were computed by Eqs. (10)-(13), where the fitting parameters λ HO , λ HH , λ OO , and λ OH were the coefficients of the interaction energies between the ith backbone dipole, � ⃗ i X (X = HN or CO), and the neighboring backbone dipole, � ⃗ j Y (Y = HN or CO) (j ≠ i), as in Eqs. (10) and (11) Here, r ij is the distance between the center of the ith backbone dipole and that of the jth dipole.
The local dipole moments were given as Kondo et al. (2019). The four parameters λ HO , λ HH , λ OO , and λ OH were determined by minimizing the square of the difference between E MM modified HB and E WH HB simultaneously for H-bond energies of 21 WH alpha values of α-helices (Kondo et al. 2019) and for those of 21 WH 3_10 values of 3 10 -helices (Kondo et al. 2022). The resulted parameters were λ HO = −11.535, λ HH = 38.517, λ OO = 16.967, and λ OO = 4.587, respectively, where the unit was all (kcal/ mol) versus E WH HB (small filled symbols) in Fig. 3 from the E WH HB values was calculated to be 0.28 kcal/mol, much smaller than the original one. It suggests that this modification of the classical MM force field could well approximate the H-bond energies provided by NFA. In particular, those for the α-helices were remarkably improved, as shown in Fig. 3. For extended structures, the amounts of changes in the atomic partial charges should be small because the dipole-dipole interaction quickly decreases in inverse proportion of the cube of the distance between the dipole pair.

Concluding remarks
We have briefly reviewed our recent computational studies of hydrogen bonds (H-bonds) in helical secondary structures of proteins, α-helix and 3 10 -helix, using a Negative Fragmentation Approach (NFA) with density functional theory (DFT). Our computation showed that the H-bond energies of the α-helix are generally weaker than those of the separated H-bonds due to the depolarized electronic structures around the carbonyl oxygen of the ith residue and the amide hydrogen of the (i + 4)th residue. The H-bond energies of the 3 10 -helix are also weaker than those of the separated H-bonds, but the effects are not so large as those for the α-helix. Whereas the adjacent C=O group is involved in the depolarization of the H-bond acceptor in α-helices, the C=O group of the H-bond pair is closer to the adjacent N-H group than to the adjacent C=O group in the 3 10 -helices. Therefore, the weak destabilization of the H-bond is attributed to the balance of the depolarization and polarization caused by the adjacent N-H group and the C=O group. Based on the findings from our computational results, a model was constructed in which the atomic partial charges of the N-H and C=O groups of the backbone peptide groups forming H-bonds are changed by the neighboring peptide groups, respectively. This modified MM model well reproduced the H-bond energies of α-helices and 3 10 -helices given by the NFA computation. We expect that this modification could lead to more reliable MD simulations in the future.