|Home | About | Journals | Submit | Contact Us | Français|
Pauling and Corey proposed a pleated-sheet configuration, now called α-sheet, as one of the protein secondary structures in addition to α-helix and β-sheet. Recently, it has been suggested that α-sheet is a common feature of amyloidogenic intermediates. We have investigated the stability of anti-parallel β-sheet and two conformations of α-sheet in solution phase using the density functional theoretical method. The peptides are modeled as two-strand Acetyl-(Ala)2-N-methylamine. Using stages of geometry optimization and single point energy calculation at B3LYP/cc-pVTZ//B3LYP/6-31G* level and including zero-point energies, thermal, and entropic contribution, we have found that β-sheet is the most stable conformation, while the α-sheet proposed by Pauling and Corey has 13.6 kcal/mol higher free energy than the β-sheet. The α-sheet that resembles the structure observed in molecular dynamics simulations of amyloidogenic proteins at low pH becomes distorted after stages of geometry optimization in solution. Whether the α-sheets with longer chains would be increasingly favorable in water relative to the increase in internal energy of the chain needs further investigation. Different from the quantum mechanics results, AMBER parm94 force field gives small difference in solution phase energy between α-sheet and β-sheet. The predicted amide I IR spectra of α-sheet shows the main band at higher frequency than β-sheet.
Pauling and Corey proposed a pleated-sheet configuration, now called α-sheet, as one of the protein secondary structures in addition to α-helix and β-sheet in 1950s.1 As shown in Fig. 1a, all of the chains in α-sheet arrange in the same way as which their main-chain carbonyl groups point in one direction and their imino groups in the opposite direction, giving rise to a strong electrostatic dipole. The inter-strand main-chain-main-chain hydrogen bonds run parallel to each other with one hydrogen bond per residue pair. The backbone dihedral angles of α-sheet are in the so called “α-extended chain” conformation, which consists of alternating right- and left-hand helical conformations, whereas both the α-helix and β-strand consist of average repeating (,ψ) angles. Yet, the occurrences of α-sheet in Protein Data Bank are not common; only seen in some special cases, i.g. K+ channel,2 or small fragments up to 4-5 residues (see Table 4 and Fig. 9 of Ref.3). Recently, α-sheet secondary structure formation was observed by molecular dynamics (MD) simulations of various amyloidogenic proteins at low pH.3-6 It was proposed that α-sheet was a common feature of amyloidogenic intermediates.6 Even at neutral pH, it was found that an engineered double mutant of human transthyretin can form α-sheet.7 However, the α-sheet structure observed from recent MD simulations is different from that proposed by Pauling and Corey in terms of hydrogen bond pattern and alignment (see Fig. 1a and b). The simulations reported a bifurcated hydrogen-bonding network (Fig. 1b). The two strands are in register in Pauling and Corey’s model, but they are shifted as observed in MD. The average backbone dihedral angles are (45°,92°) and (-87°, -49°) for αL and αR in simulation, respectively.3 The simulated α-sheet explains the hydrogen exchange data of transthyretin well.3,6
Milner-White et al. has shown the α- to β-sheet transition can occur by the flipping of alternate peptide planes.8 The peptide plane flipping in transthyretin has been described as sequential propagation like “domino effect”.6 A search within sets of closely related protein crystal structures revealed that αRαL ↔ ββ transitions occurred in 8.5% of protein families.8 The assembly of amyloidogenic yeast Sup35 assessed by scanning (atomic) force microscopy is a reminiscent of aggregation of colloidal particles in the presence of an intrinsic or induced dipole moment.9 Accordingly, it was proposed that proteins can build up a molecular dipole upon low pH conditions, resulting in attractive forces for linear self-assembly during amyloidosis.6,9 The molecular dipole model is supported by the fact that amyloid fibrils exhibit a strong tendency for self orientation and seem to form more efficiently in a magnetic field.10,11
In this work we use quantum mechanical calculations, with the aim to avoid possible force field bias in commonly used simulation packages, to explore the stability of different α-sheets in solution phase. The α-sheets we built consist of two capped dialanine peptide strands. We want to address the following questions: will the α-pleated sheet proposed by Pauling and Corey1 be stable after full optimization? How is its stability relative to MD-observed α-sheet conformation with a different hydrogen bond pattern? What is the free energy of each system? We also predict the amide I infrared (IR) spectra of the final minimized structure in solution. Our IR analysis is different from the traditional way which usually takes longer chains and multiple strands without full geometry optimization.12 This traditional approach is suitable to study regular secondary structures such as β-sheet and α-helix to investigate the structure-spectrum relationship. For the untypical conformations such as α-sheet without any available experimental spectral data, our aim is to provide clues for seeking appropriate detecting method. The recent work of Torii13 on the spectral feature of some untypical conformations used the peptide chain GNNQQNY conformation in microcrystal14 which is closely related to amyloid fibril and a fragment of Aβ(1-42) fibril.15 Unlike our solution phase study, the structural models in Torii’s work are in fact the α-sheet conformations in fibril state, which may or may not be good representatives of those in solution phase when only constrained optimizations were performed.
The sequence of a single strand is Acetyl-(Ala)2-N-methylamine. Three models of two-strand sheet were built using Maestro,16 totally sixty-four atoms in each model system. Model I was built into a standard anti-parallel β-sheet using transthyretin (PDB entry: 1TTA17) residues 107-110 and residues 117-120 as template for backbone atoms. Model II (Fig. 1a) and Model III (Fig. 1b) were built as parallel α-sheets based on Pauling and Corey’s model1 and MD simulations,3,4,7 respectively. At the building stage, the two strands were initially arranged in register. For Model I, the backbone dihedral angles were set to (-123°, 136°).8 The GB/SA implicit solvation model18 and Amber parm94 force field19 incorporated in MacroModel20 were employed. The backbone torsion angles were restrained with the force constant equal to 9999.0 kJ/(mol radius2). We carried out energy minimization using MacroModel. The convergence criterion of conjugate gradient was either the root-mean-square gradient ≤ 0.001 kJ/(mol Å) or the maximum iterations of 5000, whichever came first. For Model II and Model III, their backbone dihedral angles were set to alternate αR (-87°, -49°) and αL (45°, 92°), which are the average values observed in molecular dynamics.3 In addition to the above minimization protocol, the distance between H of backbone NH and O of backbone CO was restrained to 2.1 Å ± 0.2 Å (tolerance) with a force constant of 100.0 kJ/(mol Å2) to obtain the bifurcated hydrogen bond pattern in Model III, whereas the hydrogen bond (in green of Fig. 1b) at the center of the sheet was under more stringent restraint (2.1 Å ± 0.1 Å) and its force constant was set to 1000.0 kJ/(mol Å2). It is worth to note that the desired hydrogen bonding pattern was automatically obtained using the backbone torsion angle restraint alone for Model I and II, however, both the distance and torsion angle restraints were needed for Model III.
After the building stage, the models were subject to quantum mechanical (QM) geometry optimization. Table I summarizes the geometry optimization procedures for each model system: three stages of geometry optimization for Model I and four stages for Model II & III. The geometry optimizations were carried out using the density functional theory (DFT) method and the B3LYP exchange and correlation functionals21-23 with 6-31G* basis set in Jaguar.24 The accuracy level of SCF (self-consistent-field) calculation was set to ultrafine and the default convergence criteria were used. The solvated molecular system was treated with a self-consistent reaction field method, using Poisson-Boltzmann solver25,26 implemented in Jaguar. The dielectric constant of water was set to 80.37. The probe radius was 1.4 Å. For Model II and III, the backbone dihedral angles and the hydrogen bond distances, d(N)H...(C)O = 2.1 Å, were constrained both in vacuum and in water in stage 1 and 2 to ensure that the geometry of backbone is well preserved when the side chains are minimized. The distance constraint was released in stage 3, and finally the geometry was optimized in water without any constraint. The protocol for Model I was similar to that of Model II & III, but the distance constraint was not applied at all. Single-point energy calculation of the final optimized structure was done at the B3LYP level in solution phase with cc-pVTZ basis set. Zero-point energies and the system’s thermochemical properties were calculated at the B3LYP/6-31G* level. The frequency scale factor was 0.9614.27 The calculations of infrared (IR) intensity and frequency were also at the B3LYP/6-31G* level based on the optimized structure.
To compare with AMBER parm9419 force field, we solvated each built models in an explicit water (TIP3P)28 box with a distance of 12 Å between the wall of the box and the closest atom in the solute. Similar to the protocols of quantum mechanical geometry optimization, we employed stages of minimization. First, the conformations of water molecules were minimized for 1000 steps (50 steps of steepest descent and 950 steps of conjugate gradient) when the conformation of solute was restrained (force constant = 9999 kcal/(mol Å2)). Second, the whole solvated system was minimized for 1000 steps when the main-chain dihedral angles of solute were restrained (the relative weighting of torsion energy term = 1000). At last, when all of the restraints were released, the solvated system was minimized again for another 1000 steps. The nonbonded interactions were truncated with a 15 Å cutoff during each stage of minimization. Then the water molecules were removed in the following energy and solvation free energy calculations. The gas phase energy of the protein was calculated with a 999 Å cutoff of long-range nonbonded interaction. Delphi29 was employed to calculate the polar solvation free energy and MSMS30 was used to obtain the solvent accessible surface area (SASA). The probe radius was 1.4 Å in both Delphi and MSMS. For the Delphi calculation, the scale (the reciprocal of the grid spacing) was 2 grids/Å, the molecule filled 80% of the grid box, 1000 iterations were performed to ensure the maximum change in potential was less than 0.001 kT/e, and the internal and external dielectric constant was set to 1.0 and 80.0, respectively. The PARSE31 van der Waals radii and AMBER parm94 charges were used. The nonpolar solvation free energy was estimated by γ SASA + b, where γ and b are 0.00542 kcal/mol · Å2 and 0.92 kcal/mol, respectively.
To describe the overall twisting of sheet, we first defined a series of backbone direction vectors which start from the middle point of the first peptide bond and end at the middle point of the second peptide bond and so on. We defined the twist angle between pairs of adjacent strands in the sheet as the acute angle formed by the corresponding backbone direction vectors of two strands. The average twist angle within a sheet indicates the overall twisting of the sheet.32
The structure of β-sheet (Model I) does not change significantly during all of the stages of optimization (Fig. 2). The initially presented hydrogen bonds are well preserved. In Table II, we summarize the all-atom root mean square deviation (rmsd) between the optimized structure in each minimization stage and the structure of its previous stage as well as the hydrogen bond distance ((N)H...(C)O) in the final optimized structure. For the β-sheet, the all-atom rmsd between the structure (in Fig. 2a) resulted from the optimization in vacuum with backbone dihedral angles constrained and the built structure is 0.15 Å. In stage two and stage three, the rmsds of the optimized structure with respect to the one in the former stages are 0.14 Å and 0.22 Å, respectively.
The hydrogen bond pattern of Pauling and Corey’s α-pleated sheet (Model II) is also well maintained after four stages of minimization (Fig. 3). In solution phase the constraints were gradually removed, and the last step of minimization was without any constraint. The all-atom rmsd between the conformation (in Fig. 3a) obtained from the optimization in vacuum with both backbone dihedral angles and hydrogen bond distances constrained and the built structure is 0.39 Å, which is mainly due to the rotation of hydrogen atoms of the methyl groups of the side chains. In the second, third, and forth stage, the all-atom rmsds of the optimized structure relative to the structure of its previous stages are 0.15 Å, 0.17 Å, and 0.55 Å, respectively.
Although Model II and III are both built into α-sheets with the same backbone dihedral angles and have the same structural characteristics of which the main-chain -NH groups align on one side of the sheet and the -CO groups on the other side, their hydrogen-bonding networks are completely different. Model III has the bifurcated hydrogen bond pattern according to the molecular dynamics simulation of Armen et al.,3,6 therefore, there are five hydrogen bonds (Fig. 1b) in contrast to three hydrogen bonds in Model I (Fig. 2a) and Model II (Fig. 1a). Moreover, unlike Model II of which the two chains are in register, in Model III one strand has to shift by half a residue relative to the other strand in order to form the bifurcated hydrogen bonds, leading to the main-chain -CO groups to be located in the middle of two -NH groups of the other strand. Such shift may cause steric hindrance though Model III exhibits more hydrogen bonds than Model II. As seen in geometry optimization, the bifurcated hydrogen-bonding network is unendurable. When the (N)H···(C)O distance constraints were removed in the solution phase in stage three, one pair of bifurcated hydrogen bonds was lost (compare Fig. 4c with Fig. 4a and 4b). When the backbone dihedral angle constraints were also released in the final stage of optimization, the whole structure was deformed with only two hydrogen bonds left as shown in Fig. 4d. As pointed out by Scheiner,33 the inter-strand main-chain hydrogen bonds in extended sheets are weaker than helices and their strength is very sensitive to the backbone dihedral angles, even when the hydrogen bond geometry itself does not change from one conformation to another.
The backbone dihedral angles of the final optimized structure are presented in Table III. The backbone dihedral angles of the final optimized β-sheet structure only slightly deviate from the initial values of the built structure with the unsigned average difference of 5.4°. The largest deviation comes from 1 of the second chain, changed from -123.0° to -135.3°, but it is still within β-sheet region. Pauling and Corey described three α-sheets:1 flat, 7°-rotation, and 20° with dihedral angles of (R, ψR, L, ψL) = (-73°, -82°, 73°, 82°), (-69°, -74°, 60°, 83°), and (-97°, -66°, 56°, 106°), respectively. The average α-sheet (, ψ) angles observed from the molecular dynamics simulations by Armen et al. is αR(-87°±7°, -49°±4°) and αL(45°±8°, 92°±28°),3 which is our starting point to build the α-sheet. After the quantum mechanical geometry optimization, we have found that three out of four αL (, ψ) angles of Model II and III are in the range of (~50°, ~40°-~50°) as shown in Table III. Recent studies on alanine dipeptide conformation in water at the MP2/6-31G** level have reported αL = (59.4°, 41.4°),34 which is in good agreement with ours. This range of αL is consistent with experimentally observed bifurcation of α-sheet hydrogen bonding illustrated by Armen et al. (Fig. 9 of Ref.3). For example, the αL of Asn 608 in the NMR structure of janus-faced C2B domain of rabphilin (PDB entry: 3RPH) is (59.9°, 45.4°); the αL of Asn 207 in the crystal structure of C2 domain of snaptotagmin I (PDB entry: 1RSY) is (45.5°, 56.1°); and the αL of Asn 72 in the crystal structure of calcium phospholipid binding domain of cytosolic phospholipase A2 (PDB entry: 1RLW) is (53.2°, 52.4°). Our optimized αL dihedral angles are also close to those observed in the crystal structure of a capped tripeptide with features of α-pleated sheet, in which the dihedral angles of αL are ( = 67.3±5°, ψ = 44.1±5°).35 The αR dihedral angles of this tripeptide is ( = -75.1±5°, ψ = -25.8±6°).35 The MP2/6-31G** level geometry optimization of alanine dipeptide conformation in water gave αR = (-70.5°, -32.1°).34 The α-helical region defined in PROCHECK36 is ( = -65.3±11.9°, ψ = -39.4±11.3°). As seen in Table III, two out of four αR dihedral angles of our optimized Model II and III are in line with literature, while the other two deviate. It appears (in Fig. 3d) that even though the backbone dihedral angles of chain 1 in Model II are not within the literature reported αR and αL region, their influence on the orientation of peptide plane is small. Consequently, all of the backbone -CO groups still point to the left and -NH groups towards the right, resulting in a well maintained α-sheet structure. However, it seems that the second peptide plane of chain 1 (the left strand) in Model III (in Fig. 4d) tends to flip back to β-sheet orientation when ψ1 = -13.1° (in Table III), leading to the distortion of the chain and the disruption of the bifurcated hydrogen bond pattern.
Table IV lists the relative energy with respect to the optimized β-sheet conformation. In gas phase, β-sheet is favored by ca. 25-30 kcal/mol, compared with the α-sheets (Model II and III), while the energy difference between the β-sheet and the α-sheets in solution phase drops to only ca. 12 kcal/mol because water plays an important role in stabilizing the α-sheets by providing more favorable solvation energy contributions: 54.73 kcal/mol and 49.97 kcal/mol for the α-sheets, compared with 36.16 kcal/mol for the β-sheet. The dipole moments (Table IV) of the α-sheets are 1.4-1.7 times larger than that of the β-sheet, which is expected because the α-sheet backbone -CO groups are aligned in one direction with the -NH groups in the opposite direction. As a result, the α-sheets are more favorable in polar solvents, such as water, than in gas phase. Since the backbone -CO groups and -NH groups are better aligned in Model II final optimized structure (Fig. 3d) than that of Model III (Fig. 4d), the dipole moment of Model II is larger than Model III and its solvation energy is more favorable as well, whereas Model III is more favored than Model II in gas phase because Model II has larger nuclear-nuclear repulsion and more unfavorable two-electron term than Model III (See Supplementary Material for each component of the total energy).
We also calculated the free energy for each model after the final stage of optimization at 1.0 atm and 298.15 K in solution, including the zero-point energies, thermal, and entropic contributions (Table V). Model II has the lowest entropy (in Table V) among the three models due to the alignment of -CO groups and -NH groups. The low entropy and less favorable solution phase energy of Model II leads to an elevated free energy. Model III, though started from an α-pleated sheet, was distorted after the unconstrained optimization, thereby it has a larger entropy than Model II. The overall conclusion does not change from relative energy to relative free energy: the model β-sheet is more stable than the model α-sheets by ca. 12-13 kcal/mol in solution phase. When the chain goes longer, the favorable solvation energy relative to the increase in internal energy of the chain will change. Further calculations of longer chains and multiple strands are needed to test the possibility of the conformational conversion between β-sheet and α-sheet.
To compare with common molecular mechanics (MM) force field, we chose AMBER parm94 as an example because we reported simulation results of amyloidogenic proteins using this force field.7 We have observed marked discrepancies between the MM and DFT minimized structures. As shown in Fig. 5, both Model II and Model III have tendencies to convert to β-sheet via peptide plane flipping, especially the C-terminal peptide plane of chain 2 (the right chain) of Model II and the C-terminal peptide bonds of both chains of Model III, while the β-sheet conformation itself is well preserved. The two chains of Model II appear to slightly shift relative to each other, deviating from the aligned conformation as seen in DFT minimized structure (Fig. 3d). Accompanying these changes, the pattern of hydrogen bonding is also altered. As a result, the conformations of Model II and Model III become alike; the all-atom RMSD between Model II (Fig. 5b) and Model III (Fig. 5c) after MM minimization is only 0.78 Å, compared with 1.20 Å for DFT. The N-terminal bifurcated hydrogen bonds of Model III are kept while others disappear because of the flipping of peptide plane. We present the relative molecular mechanics energy and solvation free energy in Table VI. The energy differences between α-sheet and β-sheet in gas phase calculated with AMBER parm94 are significantly less than those calculated using quantum mechanics, down from ~25-30 kcal/mol to ~15-18 kcal/mol. The solvation free energies given by Delphi and estimated from SASA are close to the QM results. The resulting solution phase energy difference between α-sheet and β-sheet is only 1-4 kcal/mol by MM force filed, compared with ~12 kcal/mol for DFT.
The ultimate way to study the conformational change between β-sheet and α-sheet using computational approach is to carry out ab initio molecular dynamics simulations of amyloidogenic fragments in explicit water until the simulation reaches statistical convergence. Yet, it is beyond currently available computer power. Our present work should be taken as preliminary results. The distorted α-sheet of Model III in the final optimization stage could be due to the short length of our model. The shift of one dipeptide strand relative to the other in Model III α-sheet may cause destabilization. In a larger sheet, this kind shift may be better accommodated.3,6,7
Amide I IR spectral features are commonly used for determination of protein secondary structure along with other data such as circular dichroism. The assignment of amide I position of native proteins containing β-sheet is around 1633 cm-1.37 Aggregated proteins often show a maximum of IR intensity near 1620 cm-1,37-39 a characteristic of extended β-sheets with smaller twist angles than native proteins.39 It is interesting to predict the amide I band of our model systems. Our aim is not to analyze the spectral feature of MD-generated structures because those structures are likely to be biased by the particular force field used by the authors. Rather, we investigate the IR spectra of the quantum mechanical optimized structure in solution phase without any empirical parameter. For the untypical structure in this study, we believe that IR analysis based on the fully optimized structure in solution phase is the best way to truly reflect the experimental condition without any bias given the present computer power. Therefore, the results are meaningful for seeking appropriate detecting method for α-sheet. Since Model III became distorted after geometry optimization, we restricted the amide I band calculations to Model I and Model II. The solvent effect was taken into account by the continuum solvation model.25,26 As shown in Fig. 6a, the main band of anti-parallel β-sheet (Model I) is at 1620 cm-1, the second peak is at 1634 cm-1, and the third maximum is at 1655 cm-1. Besides the sensitivity on the secondary structure, the positions of amide I bands have been found to depend on solvent, number of strand per sheet, and the twisting of sheet.12,40 The overall twisting of the sheet is 9.1°, which is substantially less twisted than native proteins (30° - 40° summarized in Table 1 of Ref.39). The predicted amide I spectrum of our two-strand anti-parallel β-sheet exhibits mixed features of β-sheet in native proteins (1634 cm-1) and in aggregated states (1620 cm-1). The main band of our Model II α-sheet shifts to 1653 cm-1 (Fig. 6b). The side bands at 1629 cm-1 and 1681 cm-1 show similar intensities to each other. Following the same definition of twisting for β-sheet, the overall twisting for Model II α-sheet is 8.7°. Two major factors of α-sheet can be attributed to the frequency shift, strong permanent dipole (Table IV), and the “repulsive” electrostatic environment due to the parallel side-by-side backbone C=O alignment.
Torii, who modeled the anti-parallel β-sheet as extended state and the α-sheet as (αL=60° to 80°, αR= -60° to -80°), also predicted that the amide I frequency shifted from below 1620 cm-1 of β-sheet to above 1670 cm-1 of α-sheet (in 2H2O) and attributed the shift to the dependence of the diagonal force constants of individual peptide groups and off-diagonal vibrational coupling constants between the covalently-bonded peptide groups on the (,ψ).13 Although Torii’s calculations and ours are different in backbone dihedral angles, basis sets, the treatment of solvent effect, and most importantly, Torii took the α-sheet conformation in fibril state and treated it as if it is in solution, the consensus is that α-sheet amide I band locates in higher frequency region than that of β-sheet.
We have investigated the stability of β-sheet and two conformations of α-sheet in solution phase using quantum-mechanical approach. We have found that β-sheet is the most stable conformation, while the α-sheet proposed by Pauling and Corey has 13.6 kcal/mol higher free energy than the β-sheet. The α-sheet that resembles the structure observed in MD simulations of amyloidogenic proteins at low pH becomes distorted after stages of geometry optimization in solution. Whether the α-sheets with longer chains would be increasingly favorable in water relative to the increase in internal energy of the chain needs further investigation. Different from the DFT results, AMBER parm94 force field gives small difference in solution phase energy between α-sheet and β-sheet. The predicted amide I IR spectra of α-sheet shows the main band at higher frequency than β-sheet.
This work is partially supported by NIH (1R15 AG025023-01 to S. H) and National Science Foundation Major Research Instrumentation (DBI-0320875). The molecular graphics images were generated using the Chimera41 package from the Computer Graphics Laboratory, University of California, San Francisco, CA (supported by NIH P41 RR-01081). We thank the Scientific Computing and Visualization group at Boston University for providing part of the computational resources. Special thanks to the support team of Schrödinger, LLC for their quick and detailed response to our questions. S. Huo is grateful to Dr. Hajime Torii for helpful discussions.