Density functional study of molecular interactions in secondary structures of proteins.

Proteins play diverse and vital roles in biology, which are dominated by their three-dimensional structures. The three-dimensional structure of a protein determines its functions and chemical properties. Protein secondary structures, including α-helices and β-sheets, are key components of the protein architecture. Molecular interactions, in particular hydrogen bonds, play significant roles in the formation of protein secondary structures. Precise and quantitative estimations of these interactions are required to understand the principles underlying the formation of three-dimensional protein structures. In the present study, we have investigated the molecular interactions in α-helices and β-sheets, using ab initio wave function-based methods, the Hartree-Fock method (HF) and the second-order Møller-Plesset perturbation theory (MP2), density functional theory, and molecular mechanics. The characteristic interactions essential for forming the secondary structures are discussed quantitatively.

Proteins are essential macromolecules in biology [1], and their functions depend on their three-dimensional structures. Secondary structures, such as α-helices and β-sheets, are important components of the protein architecture. An α-helix, a particularly rigid arrangement of the polypeptide chain, is a common secondary structure element in both fibrous and globular proteins [2]. The hydrogen bonds of an α-helix are arranged such that the peptide N–H bond of the nth residue points along the helix toward the peptide C=O group of the (n−4)th residue. This results in a strong hydrogen bond that has the nearly optimum N...O distance of 2.8 Å. A β-sheet is also a common structural motif in proteins [3]. It utilizes the full hydrogen bonding capacity of the polypeptide backbone to form a rigid structure. In the β-sheet, however, hydrogen bonding occurs between protein strands, rather than within a strand as in α-helices. β-Sheets are classified into two types: antiparallel β-sheet, in which neighboring hydrogen bonded strands run in opposite directions, and parallel β-sheet, in which the hydrogen bonded strands are aligned in the same direction. The conformations in optimally hydrogen-bonded β-strands are different from those in a fully extended polypeptide, and have a pleated edge-on appearance, called a “pleated sheet”.

Molecular interactions, in particular hydrogen bonds, play significant roles in the formation of secondary structures. Precise and quantitative estimations of these interactions are required to understand the principles underlying the formation of three-dimensional protein structures.

Theoretical studies have been conducted on the hydrogen bonding interactions in secondary structures, such as α-helices and β-sheets, at different levels of theory. Yoda et al. compared the secondary structure characteristics of commonly used force fields, using molecular dynamics simulations with explicit water molecules. They showed that the molecular dynamics simulations with AMBER ff94 [4] and ff99 [5] were in remarkable agreement with experimental data for α-helical polypeptides, but provided α-helices for β-hairpin polypeptides [6,7]. Wieczorek and Dannenberg investigated the H-bond cooperativity and energetics of α-helices, and suggested that various factors contribute to their stability [8-10]. Morozov et al. assessed the origin of cooperativity in the formation of α-helices [11]. Zhao and Wu performed theoretical calculations on α-helix and β-sheet models, constructed by using a simple repeating unit approach method, to investigate the cooperativity in the formation of α-helices and β-sheets [12,13]. Their computations suggested that the formation of hydrogen bond networks in β-sheets does not occur with significant cooperativity. Parthasarathi et al. studied the hydrogen bonding interactions in α-helix and β-sheet models, using the atoms-in-molecule approach [14].

In the present study, we have investigated the molecular interactions in α-helices and β-sheets, using molecular orbital methods (the Hartree–Fock (HF) method, the second-order Møller-Plesset perturbation theory (MP2), and density functional theory (DFT)) and molecular mechanics. The characteristic interactions essential for forming the secondary structures are discussed quantitatively.

Materials and Methods

Model construction

The models of a β-strand and an α-helix were constructed by using poly-alanine amino acids capped with an acetyl group (ACE) and an N-methyl amide group (NME), denoted as ACE–(Ala)–NME (n = 1–5 for the β-strands and n = 3–8 for the α-helices), as shown in Figure 1. The parallel and antiparallel β-sheet models were constructed from a dimer of β-strands aligned in parallel and antiparallel manners, respectively. The β-sheet models have 2–6 hydrogen bonds, while the α-helix models have 1–6 hydrogen bonds. The φ and ψ angles of the parallel and antiparallel β-sheet were set to their ideal values, as described in biochemistry textbooks [15,16] (φ = −119.0° and ψ = +113.0° for the parallel β-sheet models and φ = −139.0° and ψ = +135.0° for the antiparallel β-sheet models), because β-strands usually form a pleated sheet, although the structures of the β-strands were constrained to be planar in the previous studies [12,14]. The backbone torsion angles for the α-helix models were set to φ = −57.0° and ψ = −47.0°.

Figure 1

Models of antiparallel and parallel β-sheets: (ACE-(ALA)-NME)2 and α-helix: ACE-(ALA)-NME.

Theoretical calculations

All calculations were performed on the models of the secondary structures with the Gaussian09 program packages [17]. To determine an appropriate theoretical method for the description of the electronic structures of ACE–(Ala)–NME, we assessed the validity of eleven DFT exchange–correlation functionals: B3LYP (Becke three-parameter hybrid functional combined with Lee–Yang–Parr correlation functional) [18], BLYP (Becke exchange functional combined with Lee–Yang–Parr correlation functional) [19,20], PW91 (Perdew and Wang’s 1991 functional) [21], PBE0 (hybrid funtional by Adamo with 1996 pure functional of Perdew, Burke, and Ernzerhof) [22], M06 (Minnesota 2006 functional) [23], M06-2X (Minnesota 2006 functional with 54% Hartree-Fock exchange) [23], LC-ωPBE (long range-corrected version of Perdew, Burke, and Ernzerhof functional) [24-26], CAM-B3LYP (Handy and coworkers’ long range corrected version of B3LYP using the Coulomb-attenuating method) [27], B971 (Handy, Tozer and coworkers modification to Becke’s 1997 functional) [28], B98 (Becke’s 1998 functional) [29], and B97D (B97 functional with Grimme’s D2 dispersion schemes) [30]. The computational results were compared to those at the MP2 level of theory [31]. We also compared them to those estimated by the classical force field: the AMBER force field (ff99SB [5,32]). The 6-31+G(d) basis sets were used for all calculations. The peptide models were optimized with the constraint of the φ and ψ torsion angles at the B97D/6-31+G(d) level of theory (We also performed the geometry optimization of the peptide models at the CAM-B3LYP/6-31+G(d) level of theory, and the results were similar to each other: the B97D and CAM-B3LYP functionals showed stabilization energies (SEs) close to those of the MP2 method, among the DFT exchange–correlation functionals, even when the models were optimized at the CAM-B3LYP/6-31+G(d) level of theory, as shown in Supplementary Tables S1 and S2). In the present study, we focused on hydrogen bonding energies of α-helices and β-sheets in the gas phase. Solvation effects are important because proteins are solvated in the aqueous solution. However, evaluation of solvation effects is hard and time-consuming by QM simulations because the dynamics of water molecules should be considered. It is also required to understand the intrinsic properties forming secondary structures with unsolvated models before the effects of solvation are evaluated. In addition, hydrogen bonding energies in the gas phase are useful for the improvement of the force fields.

Evaluation of stabilization energies of the β-sheet models

The stabilization energies (SEs) obtained from the formation of the parallel and antiparallel β-sheet models were evaluated with the following equation:

where Eβ-sheet is the total energy of the β-sheet model, and Epeptide1 and Epeptide2 are the total energies of the individual energies of the interacting peptides consisting of the β-sheet model. The SE values were corrected for the basis set superposition error (BSSE) by the counterpoise method of Boys and Bernardi [33].

Evaluation of relative energies of the α-helix models

The SEs for the formation of an α-helix cannot be obtained directly, because the hydrogen bonding donor and acceptor are covalently linked. Here, to estimate the stabilization energies for the formation of the α-helix model, we computed the relative energies (REs) as the difference in the total energies between the α-helix model, Eα-helix, and the same sequence in the extended conformation, Eextended, with the following equation:

Results and Discussion

Validity of DFT exchange–correlation functionals for the evaluation of hydrogen bonding interaction energies of secondary structures

Post HF calculations are required for the highly accurate evaluation of weak interactions, such as hydrogen bond energies; however, these calculations are difficult to perform, due to expensive computational costs. In recent years, DFT has become accepted as an alternative approach for the post HF methods. In previous studies [34,35], we showed the importance of assessing the validity of various DFT exchange–correlation functionals. To determine an appropriate DFT exchange–correlation functional, we computed the SEs for the formation of the parallel and antiparallel β-sheet models with 11 exchange–correlation functionals, the HF method, the post HF (MP2) method, and the classical AMBER ff99SB force field. The computed SEs for the antiparallel and parallel β-sheet models are summarized in Tables 1 and 2, and the increments of the computed SEs are listed in Tables 3 and 4, respectively. The DFT calculations provided remarkably different SEs, in the range of about 5 kcal mol−1 for the alanine monomer, ACE-(Ala)-NME, and in the range of 20 kcal mol−1 or larger for the alanine pentamer, ACE-(Ala)5-NME.

Table 1

Stabilization energies (SEs) of antiparallel β-sheet modelsa,b

Method	Basis sets	Number of hydrogen bondsc
		2	3	4	5	6
HF	6-31+G(d)	14.7	16.3	23.8	25.5	33.7
MP2	6-31+G(d)	18.7	21.9	32.6	35.9	46.8
B3LYP	6-31+G(d)	16.1	17.8	26.0	28.0	36.7
B971	6-31+G(d)	17.8	20.2	29.4	32.1	41.6
B97D	6-31+G(d)	19.5	22.6	33.6	37.1	47.8
B97D	6-31++G(d,p)	19.9	23.3	34.6	38.3	49.3
B98	6-31+G(d)	17.1	19.2	28.1	30.4	39.6
BLYP	6-31+G(d)	14.1	14.8	21.9	22.9	30.5
CAM-B3LYP	6-31+G(d)	18.6	21.4	31.2	34.2	44.3
LC-ωPBE	6-31+G(d)	16.6	18.4	27.3	29.3	38.6
M06	6-31+G(d)	20.0	22.9	34.7	37.4	49.1
M06-2X	6-31+G(d)	20.8	24.2	36.2	39.5	51.6
PBE0	6-31+G(d)	17.8	20.1	29.4	32.1	41.7
AMBER99SB		19.1	24.1	35.3	40.1	51.3

Stabilization energies are listed in kcal mol−1.

Bold fonts mean the results computed with the MP2 method and the DFT methods that provided the results close to those by the MP2 method.

The number of hydrogen bonds is calculated by using the number of alanine residues, n, through the equation: Number of hydrogen bonds = n +1.

Table 2

Stabilization energies (SEs) of parallel β-sheet modelsa,b

Method	Basis sets	Number of hydrogen bondsc
		2	3	4	5	6
HF	6-31+G(d)	10.5	16.0	20.0	25.2	29.2
MP2	6-31+G(d)	13.9	21.3	28.1	35.3	42.1
B3LYP	6-31+G(d)	11.2	17.0	21.6	27	31.6
B971	6-31+G(d)	12.9	19.4	25.0	31.1	36.7
B97D	6-31+G(d)	14.5	21.6	29.0	35.8	43.0
B97D	6-31++G(d,p)	14.9	22.4	30.0	37.2	44.6
B98	6-31+G(d)	12.2	18.4	23.6	29.5	34.6
BLYP	6-31+G(d)	9.2	13.9	17.3	21.7	25.1
CAM-B3LYP	6-31+G(d)	13.6	20.6	26.8	33.4	39.4
LC-ωPBE	6-31+G(d)	11.7	17.8	22.9	28.6	33.7
M06	6-31+G(d)	14.9	22.5	30.2	37.5	45.0
M06-2X	6-31+G(d)	15.7	23.6	31.5	39.1	46.8
PBE0	6-31+G(d)	12.8	19.3	24.9	31.0	36.6
AMBER99SB		15.8	24.5	32.9	41.1	49.5

Stabilization energies are listed in kcal mol−1.

Bold fonts mean the results computed with the MP2 method and the DFT methods that provided the results close to those by the MP2 method.

The number of hydrogen bonds is calculated by using the number of alanine residues, n, through the equation: Number of hydrogen bonds = n +1.

Table 3

Increments of stabilization energies (SEs) in antiparallel β-sheet modelsa,b

Method	Basis sets	Number of hydrogen bondsc
		2 -> 3	3 -> 4	4 -> 5	5 -> 6
HF	6-31+G(d)	1.6	7.5	1.7	8.2
MP2	6-31+G(d)	3.2	10.7	3.3	10.9
B3LYP	6-31+G(d)	1.7	8.3	2.0	8.7
B971	6-31+G(d)	2.4	9.2	2.7	9.4
B97D	6-31+G(d)	3.1	11.0	3.5	10.7
B97D	6-31++G(d,p)	3.4	11.1	3.7	11.0
B98	6-31+G(d)	2.1	8.8	2.4	9.2
BLYP	6-31+G(d)	0.7	7.1	1.0	7.6
CAM-B3LYP	6-31+G(d)	2.8	9.8	3.1	10.1
LC-ωPBE	6-31+G(d)	1.8	8.9	2.0	9.3
M06	6-31+G(d)	2.9	11.9	2.7	11.7
M06-2X	6-31+G(d)	3.3	12.1	3.2	12.1
PBE0	6-31+G(d)	2.3	9.3	2.6	9.6
AMBER99SB		5.0	11.2	4.8	11.2

Increments are listed in kcal mol−1.

Bold fonts mean the results computed with the MP2 method and the DFT methods that provided the results close to those by the MP2 method.

The number of hydrogen bonds is calculated by using the number of alanine residues, n, through the equation: Number of hydrogen bonds = n + 1.

Table 4

Increments of stabilization energies (SEs) in parallel β-sheet modelsa,b

Method	Basis sets	Number of hydrogen bondsc
		2 -> 3	3 -> 4	4 -> 5	5 -> 6
HF	6-31+G(d)	5.6	4.0	5.2	4.0
MP2	6-31+G(d)	7.4	6.8	7.2	6.8
B3LYP	6-31+G(d)	5.8	4.6	5.4	4.6
B971	6-31+G(d)	6.5	5.6	6.1	5.6
B97D	6-31+G(d)	7.1	7.4	6.8	7.2
B97D	6-31++G(d,p)	7.4	7.7	7.1	7.5
B98	6-31+G(d)	6.2	5.2	5.8	5.2
BLYP	6-31+G(d)	4.7	3.4	4.4	3.4
CAM-B3LYP	6-31+G(d)	7.0	6.1	6.6	6.0
LC-ωPBE	6-31+G(d)	6.1	5.1	5.7	5.1
M06	6-31+G(d)	7.6	7.7	7.3	7.5
M06-2X	6-31+G(d)	8.0	7.9	7.6	7.7
PBE0	6-31+G(d)	6.5	5.6	6.1	5.5
AMBERff99SB		8.7	8.4	8.3	8.4

Increments are listed in kcal mol−1.

Bold fonts mean the results computed with the MP2 method and the DFT methods that provided the results close to those by the MP2 method.

The number of hydrogen bonds is calculated by using the number of alanine residues, n, through the equation: Number of hydrogen bonds = n + 1.

As compared with the MP2 method, the AMBER ff99SB force field and the M06-2X exchange–correlation functional overestimated the stability of the β-sheet models by about 5 kcal mol−1 for the alanine pentamer, ACE-(Ala)5-NME. However, judging from Table 3, in the antiparallel β-sheet models, the SEs computed with the M06-2X functional deviated from those obtained by the MP2 method, in that the number of hydrogen bonds changed from an odd number to an even number, while those calculated with the AMBER force field deviated, in that they changed from an even number to an odd number. The B3LYP, B971, B98, BLYP, HF, LC-ωPBE, and PBE0 calculations underestimated the stabilization energies by over 5 kcal mol−1, in comparison to the calculations by the MP2 method. In particular, it is noteworthy that the conventional B3LYP and BLYP functionals provided more than 10 kcal mol−1 lower SEs than those obtained by the MP2 method. In addition, since the increments of the SEs calculated with the B3LYP and BLYP methods were about 2–3 kcal mol−1 smaller than those computed with the MP2 methods, the discrepancies increase with the lengths of the β-strands.

The B97D and CAM-B3LYP functionals showed SEs close to those obtained with the MP2 method, among the DFT exchange–correlation functionals used in the present study. The increments of the SEs with the B97D functional were equivalent to those with the MP2 method, implying that the B97D method can be expected to provide SEs close to those computed with the MP2 method, even when the peptide lengths are elongated. However, the mean deviation of the increments with the CAM-B3LYP functional from the MP2 method was evaluated to be 0.6 kcal mol−1, and was twice as large as that obtained with the B97D method. The B97D and CAM-B3LYP functionals showed a clear improvement over generalized gradient approximation and hybrid functionals without dispersive and long-range interactions, B3LYP, B971, B98, BLYP, and PBE0. It indicates that the inclusion of van der Waals interactions is, in general, required for a quantitative description of hydrogen bonds in α-helices and β-sheets, as well as electrostatic interactions.

We also investigated the dependence of the basis set of hydrogen atoms on the SEs. As shown in Tables 1 and 2, the SEs computed with the B97D/6-31++G(d,p) method were almost the same as those obtained with the B97D/6-31+G(d) method, and the difference was within 1.6 kcal mol−1. This indicates that the effect of the polarization and diffuse functions of the hydrogen atoms is small.

In the present study, the B97D/6-31+G(d) method was chosen as the most reliable to evaluate the stabilization of the hydrogen bonding interaction energies of protein secondary structures.

Stabilization energies of the antiparallel β-sheet models

The SEs and the increments of the SEs of the antiparallel β-sheet models are listed in Figure 2. In the antiparallel β-sheet models, we found two types of stabilization, small and large stabilization. The addition of the third and fifth hydrogen bonds resulted in the small stabilization, while the large stabilization was found by the addition of the second and fourth hydrogen bonds, as illustrated in Figure 2B. This is because the odd-numbered β-sheet models formed small hydrogen bond ring structures, and the even-numbered β-sheet models formed large hydrogen bond ring structures. Although Zhao et al. [12] and Parthasarathi et al. [14] reported that large stabilization and small destabilization were found in the formation of the antiparallel β-sheet models, using the planar antiparallel sheets, our computations indicated two types of stabilization: large and small stabilization, using more realistic pleated β-sheet models.

Figure 2

Relationships of stabilization energy (A) and increment of stabilization energies (B) to the number of hydrogen bonds in the antiparallel β-sheet models. The number of hydrogen bonds is calculated by using the number of alanine residues, n, through the equation: Number of hydrogen bonds = n +1.

As compared with the B97D/6-31+G(d) calculations (QM), the AMBER force field (MM) overestimated the SEs in the odd-numbered antiparallel β-sheet models. This means that the deviation between QM and MM occurred in the evaluation of the small stabilization, as illustrated in Figure 2B, thus indicating that the overestimation of the SEs by MM is due to the underestimation of the electrostatic and van der Waals repulsion. Table 3 shows that the small stabilization computed with MM was larger than that computed by quantum chemical calculations including the DFT calculations, implying that the underestimation of the electrostatic repulsion was an intrinsic characteristic of MM. The electrostatic interaction of MM must be improved for the accurate description of the hydrogen bonding interaction energy.

Stabilization energies of the parallel β-sheet models

The SEs and the increments of the parallel β-sheet models are shown in Figure 3. The computational results using both QM and MM revealed an almost linear relationship between the SEs and the number of peptide bonds in the parallel β-sheet models, since the Pearson’s correlations were evaluated to be 0.99994 for QM and 0.99995 for MM, respectively. This is due to the formation of H-bond ring structures with almost the same size. As compared with QM, MM provided larger SEs in the parallel sheets. In addition, the increments of MM were about 1 kcal/mol larger than those of QM, indicating that the deviation of the SEs calculated with MM from QM increased with the peptide length of the β-strand.

Figure 3

Relationships of stabilization energy (A) and increment of stabilization energies (B) to the number of hydrogen bonds in the parallel β-sheet models. The number of hydrogen bonds is calculated by using the number of alanine residues, n, through the equation: Number of hydrogen bonds = n +1.

We also compared the SEs between the antiparallel and parallel β-sheet models. The antiparallel β-sheet model was more stable than the parallel β-sheet model for the same peptide length, as shown in Figures 2A and 3A. This is due to the difference in the hydrogen bond structures. The structural parameters of the hydrogen bonds in the antiparallel β-sheet models were optimized to be 1.88 Å for the O…H length, 153° for the C=O…H angle, and 165° for the O…H–N angle on average, while those in the parallel β-sheet models were 1.84 Å, 160°, and 170°, respectively. Previously, Morozov and his coworkers performed the MP2/aug-cc-pVDZ calculations on a formamide dimer, which is a model of hydrogen-bonded peptides, to investigate the orientation dependence and the distance dependence of the hydrogen bonding energies [36]. The optimized structural parameters of the formamide dimer were as follows: O…H = 1.97 Å, ∠C=O…H = 110°, and ∠O…H–N = 155°. This showed that the hydrogen bond structures of the antiparallel β-sheet models were closer to those of the optimized hydrogen-bonded formamide dimer, in comparison with the parallel models, demonstrating that the antiparallel β-sheet models had more ideal hydrogen bonds, and thus were more stable than the parallel models. This result indicates the importance of the incorporation of the anisotropic nature in hydrogen bonds, such as the C=O…H and O…H–N angles, into the force field.

Relative energies of the α-helix models

The REs and the increments of the α-helix models are shown in Figure 4. In contrast to the SEs of the β-sheet models, the increment of the interaction energies gradually increased together with the number of hydrogen bonds in the α-helices, indicating the cooperativity in the α-helix models, as illustrated in Figure 4B. The MM force field overestimated the REs of the α-helix models, as compared with the QM method. In addition, the deviations of the stability in the formation of α-helices between QM and MM were much larger than those in the formation of β-sheets. This corresponds to the fact that in folding simulations, the AMBER ff99SB force field tends to provide an α-helix for a peptide that actually forms a β-sheet in the native state [5,32,37]. It also implies that the electrostatic interaction must be improved.

Figure 4

Relationships of relative energy (A) and increment of stabilization energies (B) to the number of hydrogen bonds in the α-helix models. The number of hydrogen bonds is calculated by using the number of alanine residues, n, through the equation: Number of hydrogen bonds = n − 2.

The QM method revealed that the hydrogen bonding interaction energies of the α-helix models were less stable than those of the β-sheet models with the same number of hydrogen bonds, because the average O...H length in the α-helix models was evaluated to be 2.30 Å, and was about 0.4 Å longer than that in the β-sheet models.

Conclusion

We performed quantum chemical (QM) methods (the Hartree–Fock (HF) method, the second-order Møller-Plesset perturbation theory (MP2), and density functional theory (DFT)) and molecular mechanical (MM) calculations on models of protein secondary structures to understand how molecular interactions, in particular hydrogen bonds, function in the formation of the secondary structures. Our conclusions were as follows: (i) the comparison of the SEs of the β-sheet models to the MP2 method demonstrated that the B97D/6-31+G(d) method was the most reliable; (ii) the conventional B3LYP functional evaluated much lower SEs than the MP2 method; (iii) small and large stabilizations occur in the antiparallel β-sheet models, due to the electrostatic interactions; (iv) the AMBER force field overestimated the evaluation of the small stabilization, because the electrostatic repulsion was underestimated; (v) a linear relationship was found between the SEs and the number of peptide bonds in the parallel β-sheet models; (vi) the antiparallel β-sheet model provided a larger SE than the parallel β-sheet model for the same peptide length, because the hydrogen bond structures of the antiparallel β-sheet models were closer to those of the optimized hydrogen-bonded formamide dimer, in comparison with the parallel model; (vii) the MM force field overestimated the reaction energies of the α-helix models, as compared with the QM method, and the overestimation was consistent with the fact that the folding simulation with the AMBER ff99SB force field favored α-helices over β-sheets; and (viii) the hydrogen bonding interaction energies of the α-helix models were less stable than those of the β-sheet models, due to the longer O...H lengths in the α-helix models. The implementation of these insights will improve the force fields for accurate descriptions of protein behaviors.