Structure and Conformational Stability of Protonated Dialanine - The

Jun 24, 2010 - A systematic investigation on the structure and stability of the four conformers of the protonated dialanine cation (transA1, transA2, ...
0 downloads 6 Views 2MB Size
Download PDF
J. Phys. Chem. A 2010, 114, 7537–7543

7537

Structure and Conformational Stability of Protonated Dialanine Woon Yong Sohn† and Jae Shin Lee*,‡ DiVision of Energy Systems Research and Department of Chemistry, Ajou UniVersity, Suwon, Korea, 443-749 ReceiVed: January 14, 2010; ReVised Manuscript ReceiVed: May 7, 2010

A systematic investigation on the structure and stability of the four conformers of the protonated dialanine cation (transA1, transA2, transO1, cisA3) was performed employing the HF, MP2, and hybrid DFT methods with various basis sets which ranged from the 6-31G* to the basis set larger than the correlation consistent aug-cc-pVTZ basis set. It is found that the backbone dihedral angles and energies of the conformers are sensitive to the electron correlation level and basis set, especially manifesting slow convergence of conformer structure and energetics with basis set at the MP2 level. At the MP2 basis set limit corrected by CCSD(T) correlation effect, the lowest transA1 conformer is almost isoenergetic with the cisA3 conformer, followed by the transA2 conformer (∼0.5 kcal/mol above transA1), and, last, the transO1 conformer (∼1.2 kcal/mol above transA1). Vibrational and thermal (entropic) factors appear to have an important effect on the relative stability between conformers at room temperature, reducing the energy difference between transA1 and transA2 conformers and making cisA3 higher in energy than transA1 or transA2, which is in accord with the recent infrared multiphoton dissociation experimental data on this cation. According to the polarizable continuum model calculations, solvation of protonated dialanine in water would significantly enhance the stability of the transA2 conformer, making it most populated in aqueous solution at room temperature. Among the tested hybrid DFT methods in this study, B3LYP/6-31G* was found to be the most effective for predicting the conformational structures and relevant stability of protonated dialanine cation in the gas phase. I. Introduction Dialanine is an interesting and important model system to understand the structure and dynamics of peptides. As one of the simplest peptides having one peptide bond, it is relatively easy to explore the PES (potential energy surface) of the molecule theoretically and, in combination with experimental data analysis, understand how the peptide bond interacts with side groups (in this case the methyl group) and terminal groups to affect the overall structure of the backbone of the peptide. As a good example, one can point out recent IRMPD (infrared multiphoton dissociation) experiments on its protonated cation, protonated dialanine (Alanine-Alanine-H+), and relevant theoretical calculation studies.1,2 From the DFT (density functional theory) based quantum mechanical calculations and molecular dynamics simulation results, it was possible to understand the stability between various conformers and relevant conformational dynamics as well as spectral assignment of the vibrational modes of this cation. Indeed, the authors in previous papers were able to identify the lowest energy conformer and find that the trans form of the protonated dialanine is more stable than its cis counterpart (trans, cis refers to the orientation between carbonyl and amino group in the peptide bond). They also found that the excess proton primarily resides in terminal amino sites rather than the adjacent carbonyl group, which is able to attract the excess proton through hydrogen bonding. However, it has to be noted that many of these findings are based on theoretical calculations employing DFT methods, specifically B3LYP3-5 and BLYP3,6 methods, and considering the small energetic differences between various conformers in the protonated dialanine potential energy surface, it is not clear how structures * To whom correspondence should be addressed. † Division of Energy Systems Research. ‡ Department of Chemistry.

and stability of the protonated dialanine conformers would be affected when more accurate theoretical methods than the previously used methods (hybrid DFT methods) are employed along with large basis sets sufficient to cover the diffuse electronic motions in various conformers. Therefore, to accurately determine the structures and stabilities of various conformers of the protonated dialanine and to understand the dynamics of these conformers in the PES correctly, it appears necessary to investigate the protonated dialanine by using more sophisticated theoretical methods such as Moller-Plesset perturbation theory (MPPT)7 or coupled cluster methods8 which can account for a more accurate correlation effect on the structure and energetics of this molecular ion. In this paper we examine the equilibrium structures and conformational dynamics of protonated dialanine using the HF (Hartree-Fock), MP2 (second order Moller-Plesset perturbation theory), and several hybrid DFT methods, as well as the barrier heights between the conformers. Furthermore, we examine higher order electron correlation effects beyond MP2 on the structure and conformational dynamics of the protonated dialanine by employing the CCSD(T) (singles and doubles coupled cluster method with perturbative triples correction)8,9 method. The electronic and vibrational contributions to the stability of various conformers of the protonated dialanine will be carefully analyzed. After structures and stabilities of important (low energy) conformers of the protonated dialanine in the gas phase are determined, we investigate the structural and energetic change of these conformers in aqueous environment, which would be more realistic than the gas phase as most peptides are dissolved in water in biological systems. This paper is organized as follows. In section II computational methodologies are described in detail. In section III the results for the equilibrium structures and relative stability of the protonated dialanine conformers in the gas phase as well as in

10.1021/jp100380c  2010 American Chemical Society Published on Web 06/24/2010

7538

J. Phys. Chem. A, Vol. 114, No. 28, 2010

Sohn and Lee

Figure 1. A schematic view of protonated dialanine conformers: (top) transA1 (left) and transA2 (right); (bottom) transO1 (left) and cisA3 (right).

the aqueous solution are presented. The summary of the results of this work is given in section IV. II. Computational Details The conformer notations of the protonated dialanine used in this work, transA1, transA2, transO1, and cisA3, are the same as the conformer notations examined in previous studies.1,2 A and O distinguish the protonation site, with A representing protonation in the N-terminal amino group and O representing protonation in carbonyl oxygen in the CdO amide group, respectively. Trans and cis refer to the orientation of NH and CO groups in the peptide bond. A schematic view of the four conformers is shown in Figure 1. TransA1 and transA2 (and also cisA3) are also distinguished by the extent of rotation around the N2-C3 bond. Although these lowest energy conformers were found by extensive search of the semiempirical potential energy surface of the protonated dialanine, to check whether there are any other conformers with lower energies than these conformers, we first performed several molecular dynamics (MD) simulation calculations at 298 K using OPLS-AA force field10,11 (implemented in Gromacs program package).12,13 The configurations of MD trajectories were classified according to similar structure and eventually yielded the same conformer type corresponding to the transA1, transA2, and cisA3 conformers, although the initial structures for the MD run were chosen to belong to the same conformer type in each case. For transO1 configuration for which normal peptide bond is disrupted by proton migration, direct quantum chemical optimization from arbitrary initial structure with hydrogen attached to amide oxygen was performed without prior MD simulation (see below). Since it became evident that these four types of conformers correspond to the lowest energy conformers of the protonated dialanine, they were then optimized using three different quantum chemical methods: HF, MP2, and hybrid DFT methods which included B3LYP,3-5 BLYP,3,6 BPW91,6,14 and B3PW915,14 methods. Since basis sets can have a significant effect on the optimized structures and energetics of conformers, various bases, including 6-31G*, 6-31+G*, 6-311++G**, and correlation consistent (aug-)cc-pVXZ (X ) D, T) basis sets along with the aug′-cc-pVTZ basis set which included the diffuse functions only for the heavy atoms have been employed. To investigate the higher electron correlation effect beyond the MP2 level on the stabilities of the conformers, CCSD(T) calculations at the

MP2/aug-cc-pVTZ optimized geometries of the conformers were carried out with the cc-pVDZ, 6-31G*, 6-31+G*, and 6-311++G** basis sets. Since these calculations were carried out under frozen core approximation,17 core-correlation effect was investigated by employing the cc-pCVDZ basis set,18 which contains core-correlating functions in addition to the valenceelectron optimized cc-pVDZ set, for all electron correlation calculations. The vibrational and thermal effect on the stability of the conformers were examined at the B3LYP/6-31G* level at which simulated IR (infrared) spectra for each conformer type were also produced to compare with the experimental spectrum of this cation. To investigate the pathway of proton migration from transA1 configuration to transO1 configuration, intrinsic reaction coordinate calculation19 along with transition state optimization was performed at the B3LYP/6-31G* level. The barrier heights between the transA1 and other conformations were also evaluated for population analysis of different conformers at room temperature. Finally, to investigate the structural and energetic changes of the protonated dialanine in aqueous solution, the PCM (polarizable continuum model) method20-25 was employed, which treats water solvent as a continuous dielectric medium, at the B3LYP/6-31G* level. As in the case of the gas phase, the effect of nuclear motions on relative stability between the conformers of the protonated dialanine in aqueous solution was analyzed after vibrational frequency calculations and analysis of thermal energy contributions. All quantum chemical calculations were performed with the Gaussian03 program package.26 III. Results and Discussion The structure of a peptide can be described by the dihedral angles along the backbone of the peptide. For protonated dialanine, there are four dihedral angles along the backbone including N- and C-terminals: Ψ1 ) N1-C1-C2-N2, ω ) C1-C2-N2-C3, Φ ) C2-N2-C3-C4, and Ψ2 ) N2-C3C4-O4 (see Figure 1). In Table 1, we present these four dihedral angles of four lowest conformers of protonated dialanine optimized at various levels according to the basis set and electron correlation treatment. To evaluate how each level of theory performs in predicting the reliable geometries of this protonated peptide, we also present the average of absolute deviations |δ|av from the reference dihedral angles, which were taken as the dihedral angle values for the MP2/aug-cc-pVTZ

Structure and Stability of Protonated Dialanine

J. Phys. Chem. A, Vol. 114, No. 28, 2010 7539

TABLE 1: Differences in Dihedral Anglesa from the Reference Values According to Correlation Level and Basis Set for Protonated Dialanine Conformers transA1 Ψ1

ω

Φ

transA2 Ψ2

Ψ1

ω

HF/cc-pVDZ -10.0 -2.5 6.1 0.3 -8.3 -1.3 HF/cc-pVTZ -9.5 -2.2 6.9 -1.6 -0.3 1.6 HF/aug-cc-pVTZ -9.6 -2.0 6.9 -1.8 -6.9 2.0 BLYP/6-31G* 3.5 1.2 0.8 -0.4 -1.3 0.6 B3LYP/6-31G* 2.6 1.1 0.5 0.1 2.5 0.7 B3LYP/6-31+G* 0.7 0.7 3.5 -1.7 -3.3 3.9 B3LYP/cc-pVTZ 2.7 0.8 3.2 -1.1 3.9 4.9 B3PW91/6-31G* 3.0 1.3 0.5 0.2 -0.8 2.3 BPW91/6-31G* 3.5 1.3 0.9 -0.3 0.4 3.4 MP2/6-31G* -8.3 -1.9 0.7 0.3 -7.7 -1.6 MP2/6-31+G* -11.2 -2.3 3.0 -1.7 -8.7 -0.4 MP2/6-311++G** -4.4 -0.6 2.2 -2.0 -5.2 -1.0 MP2/cc-pVDZ -0.3 -0.2 -1.0 1.3 -1.8 -2.6 MP2/cc-pVTZ 0.0 0.1 -0.6 0.2 0.1 0.3 MP2/aug-cc-pVDZ -1.7 0.1 -0.4 -0.9 -1.6 -1.0 MP2/aug′-cc-pVTZ -0.2 0.1 0.2 0.0 -0.2 0.3 ref anglesc 173.3 176.2 197.4 173.5 168.2 161.6 a

Φ

transO1 Ψ2

Ψ1

ω

-10.2 5.4 -8.2 -1.0 -9.4 1.0 -8.4 -0.9 -9.0 0.2 -8.3 -0.8 -12.8 6.2 2.6 1.2 -8.2 6.2 2.4 1.2 -12.7 1.6 1.3 1.1 -11.0 2.5 1.2 0.8 -13.0 6.3 2.2 1.1 -15.6 6.1 2.5 1.1 -3.0 5.1 2.4 0.7 -2.1 0.7 -0.5 0.4 -0.7 1.0 0.4 0.8 -3.0 7.5 0.5 0.1 -1.5 1.2 -0.2 -0.1 1.4 -0.6 -0.6 0.3 -0.3 0.1 -0.1 0.1 296.9 171.5 176.8 178.0

cisA3

Φ

Ψ2

Ψ1

ω

5.1 6.4 6.4 -0.4 -0.5 1.7 2.0 -0.5 -0.4 0.3 1.4 1.8 -1.1 -0.3 -0.2 0.2 196.2

-0.5 -2.1 -2.2 0.0 0.4 -1.3 -0.9 0.5 0.3 -0.1 -1.7 -2.1 1.2 0.2 -0.9 -0.1 172.8

-0.8 -1.6 -2.0 -6.1 -4.7 -6.7 -6.5 -3.6 -4 1.2 0.5 1.5 1.9 0.4 1.2 0.2 133.5

-3.6 -0.4 0.1 1.2 1.4 3.5 3.6 1.4 1.2 -2.7 -1.4 -2.6 -5.6 -0.9 -1.2 -0.5 332.8

Φ

|δ|avb

Ψ2

-0.2 6.5 -1.6 3.6 -1.6 2.8 -4.1 10.2 -2.8 7.0 -4.4 6.5 -4.2 6.1 -2.6 5.9 -3.2 7.7 0.1 3.7 -0.1 0.9 1.6 -0.9 0.9 7.0 -0.4 2.5 0.0 1.5 0.2 0.4 308.7 198.5

4.4 4.0 3.2 3.2 2.6 3.4 3.3 2.8 3.3 2.5 2.3 1.8 2.3 0.9 0.9 0.2 0.0

All values in units of degree. b Average values of absolute deviations from the reference values. c MP2/aug-cc-pVTZ optimized angles.

optimized geometries of each conformer. The choice of the MP2/aug-cc-pVTZ optimized angles as the reference values was made because the dihedral angles appeared to change negligibly with basis set increase beyond the aug-cc-pVTZ set at the MP2 level. Therefore, |δ|av could be a good measure for the accuracy of a given theoretical model in determining the structure of the protonated dialanine and relevant peptides. As expected, the HF method performs most poorly among the three theoretical methods. Even with the basis set as large as the aug-cc-pVTZ, which appears to be close to the basis set limit, |δ|av amounts to more than 3° and some dihedral angles deviate from the reference values by as much as 9°. For the correlated results, although the deviation generally appears to become smaller as the basis set increases, increase of basis set does not always improve the accuracy of the optimized structures. This is well-manifested in the case of the MP2/631G* and MP2/6-31+G* results for transA1 conformer, where the dihedral angles in the latter are worse than the dihedral angles in the former. This pattern is more evident for the B3LYP/6-31G* and B3LYP/6-31+G* optimized results, where the former results are better than the latter (in terms of |δ|av) for all conformers. In fact, even at the MP2/6-31+G* or B3LYP/ 6-31+G* level, some dihedral angles are different from the reference angles by more than 10° in certain conformations. Indeed, |δ|av at the B3LYP/cc-pVTZ level (which is close to the basis set limit at the B3LYP level) is larger than |δ|av at the B3LYP/6-31G* level by almost 1° (3.3° vs 2.6°). Considering accuracy and computational feasibility, B3LYP/6-31G* (or MP2/6-31G*) appears to be a reasonable method to examine the structures of larger peptides than examined in this study. However, despite the apparent accuracy observed in the average deviation from the reference values, a caution appears necessary in assessing the accuracy of the optimized structures even at this level as deviations of dihedral angles in certain cases could be larger than 5° from the reference values, which is exemplified in the case of Φ for transA2 conformer. Finally, it is worth noting that the MP2/aug′-cc-pVTZ optimized results are almost the same as the reference (MP2/aug-cc-pVTZ optimized) results, suggesting appropriate addition of diffuse functions would be sufficient to ensure converged geometries (dihedral angles) as long as a proper number of other (polarization) basis functions are present.

TABLE 2: Relative Electronic Energies of TransA2, TransO1, and CisA3 Conformers with Respect to TransA1 Energies (in kcal/mol)b HF/cc-pVDZ HF/cc-pVTZ HF/aug-cc-pVTZ BLYP/6-31G* B3LYP/6-31G* B3LYP/6-31+G* B3LYP/cc-pVTZ B3PW91/6-31G* BPW91/6-31G* MP2/6-31G* MP2/6-31+G* MP2/6-311++G** MP2/cc-pVDZ MP2/cc-pVTZ MP2/aug-cc-pVDZ MP2/aug′-cc-pVTZ MP2/aug-cc-pVTZ MP2/CBSb

transA2

transO1

cisA3

|δ|ava

0.00 0.30 0.27 2.04 1.72 1.26 1.51 1.89 2.26 0.59 -0.24 -0.19 0.52 0.68 0.07 0.42 0.47 0.72

-0.35 -0.43 -0.35 2.70 2.34 3.18 1.21 2.16 2.31 4.92 5.59 1.82 1.51 1.97 2.73 2.46 2.42 2.24

-0.78 0.56 0.63 1.52 0.76 1.23 2.15 1.36 2.39 -1.08 -1.02 -0.81 -0.91 0.47 -0.72 0.19 0.07 0.51

1.53 1.04 1.05 0.93 0.45 0.74 1.15 0.70 1.16 1.46 1.95 0.88 0.78 0.11 0.79 0.28 0.29 0.00

a Average values of the absolute deviations of the relative energies from the MP2/CBS values for transA2, transO1, and cisA3 conformers. b Estimated MP2 basis set limit values (see the text).

In Table 2 we present the energies of each conformer at equilibrium (∆Eel) with respect to the energy of the transA1 conformer at equilibrium according to electron correlation and basis set. The values for the MP2/CBS in the last row represent the estimated basis set limit values obtained by extrapolating the correlation energies ECORR(X) with the aug-cc-pVDZ and aug-cc-pVTZ basis set by ECORR(∞) ) ECORR(X) + A/(X + 1)3 (X ) 2 for DZ, 3 for TZ)27 at the MP2/aug-cc-pVTZ optimized geometries, where ECORR(∞) corresponds to the (estimated) basis set limit correlation energies at the MP2 level. The HF contribution at the basis set limit was estimated similarly by using the X-3.4 extrapolation formula by Truhlar.28 To confirm the validity of the MP2/CBS values obtained in this way, we further performed optimizations for the transA1 and cisA3 conformer with the aug′-cc-pVTZ* basis set, which consists of the aug′-cc-pVTZ set augmented by highest polarization functions from the cc-pVQZ set (i.e., g-type functions for C, N, O and f-type functions for H). The MP2/aug′-cc-pVTZ* value of

7540

J. Phys. Chem. A, Vol. 114, No. 28, 2010

Sohn and Lee

TABLE 3: MP2 and CCSD(T) Energies of the Protonated Dialanine Conformersa transA2

b

transO1 b

basis set

MP2

CCSD(T)

∆C(T)

cc-pVDZ 6-31G* 6-31+G* 6-311++G**

0.63 0.67 -025 -0.18

0.47 0.47 -0.50 -0.36

-0.16 -0.20 -0.25 -0.18

cisA3

MP2

CCSD(T)

∆C(T)

MP2

CCSD(T)

∆C(T)

1.56 4.76 5.31 1.81

0.74 3.74 4.21 0.92

-0.82 -1.02 -1.10 -0.89

-0.86 -1.27 -1.37 -0.80

-1.20 -1.75 -1.71 -1.18

-0.34 -0.48 -0.34 -0.38

a Energies (in kcal/mol) with respect to the energies of transA1 conformer computed at the MP2/aug-cc-pVTZ optimized geometries. Differences between the MP2 and CCSD(T) relative conformer energies.

0.3 kcal/mol for cisA3 conformer energy is close to the MP2/ CBS result of 0.5 kcal/mol and implies the MP2/CBS values are reliable estimates to the basis set limit values. In Table 2 one can notice that like the variation in the case of dihedral angles in Table 1, the energy difference between conformers exhibits a wide variation according to electron correlation and basis set, especially when the basis set employed is far from convergence. It is remarkable to note that even at the MP2/6-311++G** level, cisA3 and transA2 would be predicted to be more stable than transA1, in contradiction with the MP2/CBS limit results. Likewise, the MP2/aug-cc-pVDZ optimization results manifest for cisA3 would be more stable than those for transA1 by 0.7 kcal/mol. Other than the MP2/ aug-cc-pVTZ (or MP2/aug′-cc-pVTZ) and MP2/cc-pVTZ optimizations for which |δ|av for conformer energy differences are less than 0.3 kcal/mol, B3LYP/6-31G* optimizations provide reasonably reliable results for conformer energy differences as well as geometries, ∆Eel being close to ∆Eel(MP2/CBS) within 0.5 kcal/mol on the average. Although MP2 level is often considered to be an adequate level to evaluate the correlation effect in various molecular systems including hydrogen-bonded systems,29-31 it is not clear how much further electron correlation effect beyond the MP2 level would change the conformer energies and energy differences, especially considering small energy differences between the conformers of the protonated dialanine species. For this purpose, we have computed the conformer energies at the CCSD(T) level with the cc-pVDZ, 6-31G*, 6-31+G*, and 6-311++G** basis set at the MP2/aug-cc-pVTZ optimized geometries of the conformers. The conformer energies with respect to transA1 conformer energy at the CCSD(T) level along with the corresponding MP2 results, which are summarized in Table 3, show that higher order correlation effect beyond the MP2 level could be nonnegligible, reducing the energy differences between transA1 and the other conformers (transA2, transO1, cisA3) by about 0.2-1 kcal/mol as the correlation level changes from the MP2 to CCSD(T) level. Therefore, taking into account this higher order correlation effect to the MP2/CBS values in Table 1, the two lowest energy conformers, transA1 and cisA3, appear almost isoenergetic, followed by transA2 and transO1, in order. In view of the electronic isoenergeticity of the transA1 and cisA3 conformer, the relative stability between the conformers could be determined by vibrational and thermal (entropic) effects. To explore these factors, we computed the zero point energies of each conformer and difference in thermal contributions to the energy, enthalpy, and Gibbs free energy (with respect to transA1) at 298 K employing the B3LYP/6-31G* level. The results, which are summarized in Table 4, show that the vibrational motion and entropic effect appear to play a critical role in enhancing the stability of the trans conformers over the cisA3 conformer at room temperature, although these results are based on harmonic approximations of the vibrational modes

TABLE 4: Electronic Energies (∆Eel), Zero Point Energies (ZPE), and Thermal Contributions to Energy (∆E298), Enthalpy (∆H298), and Gibbs Free Energy (∆G298) of the Protonated Dialanine Conformers at 298 Ka ∆Eelb ZPE ∆E298c ∆H298c ∆G298c

transA1

transA2

transO1

cisA3

0 0 0 0 0

0.54 0.03 0.11(0.65) 0.11(0.65) -0.13(0.41)

1.24 -0.61 -0.74(0.50) -0.74(0.50) -0.16(1.08)

0.13 0.87 0.64(0.77) 0.64(0.77) 1.93(2.06)

a Zero point energies and thermal contribution values are with respect to the values for transA1 conformer at the B3LYP/6-31G* level, in units of kcal/mol. b Eel corresponds to the MP2/CBS values in Table 2 added by CCSD(T) correction term with the 6-311++G** in Table 3. c Thermal contributions at 298 K and 1 atm. Values in parentheses represent the sum of electronic and thermal contribution.

of the protonated dialanine species. Considering electronic and thermal contributions together, transA1 is the most stable in terms of free energy, followed by transA2, transO1, and finally cisA3, in order at room temperature. The small energy difference between transA1 and transA2 forms appears to explain nonsymmetric doublet peaks near 1750 cm-1 in the IR spectrum of the protonated dialanine at 300 K observed in recent IRMPD experiment.1 As clearly seen in Figure 2 where simulated IR spectra of the transA1 and transA2 conformers (obtained at the B3LYP/6-31G* level with a regular scaling factor of 0.961432 and fwhm of 4 cm-1) are presented and compared to the experimental spectrum of the protonated dialanine, one of the most distinctive IR features between the transA1 and transA2 conformer appears to be the degree of asymmetry in the intensities of the doublet peaks near 1750 cm-1 which are related to the motions of the two carbonyl groups in the protonated dialanine cation. The exhibition of clear nonsymmetric doublet peaks near 1750 cm-1 in the IRMPD spectrum of the protonated dialanine appears to suggest significant population of transA2 conformer should be present at room temperature, which is in accord with our high level ab initio calculation results, though the peak intensities in the IRMPD spectrum cannot be directly correlated to the population of the existing conformers due to the internal vibrational redistribution in the IRMPD process.1 One important subject to be dealt with concerning conformational dynamics of protonated dialanine is to determine accurate barrier heights for conformational changes between various conformers, including the possible migration of proton from the terminal amine group to the adjacent oxygen in the amide carbonyl group by making the N-H-O hydrogen bond. The possibility of proton migration is an important question in view of certain fragmentation patterns in the mass spectrometry experiments, which suggested a strong possibility of relatively free migration of proton to other sites from basic sites.33-38 For this purpose, transition state optimizations were performed for

Structure and Stability of Protonated Dialanine

J. Phys. Chem. A, Vol. 114, No. 28, 2010 7541 TABLE 5: Barrier Heightsa for Conformational Changes between the Protonated Dialanine Conformers reactant transA1 B3LYP/6-31G* MP2/6-311++G**

product transA2 1.78(0.96) 0.87

cisA3 b

28.55(26.49) 25.32

transO1 b

3.63 2.90

a Electronic energies (in kcal/mol) of the transition state with respect to the ground electronic energies of the transA1 conformer. b Barrier heights for conformational change in aqueous solution (see the text).

Figure 2. IR spectrum of protonated dialanine in the gas phase: (a) simulated IR spectra of the transA1 and transA2 conformer; (b) experimental IRMPD spectrum of protonated dialanine from ref 1. Reproduced by permission of the PCCP Owner Societies 2004.

transA1-transA2, transA1-cisA3, and transA1-transO1 conformational changes at two different theoretical levels: B3LYP/ 6-31G* and MP2/6-311++G** level. The calculations at the latter theoretical level were performed to assess the dependence of the barrier heights on the theoretical model employed. In Table 5 barrier heights for conformational changes corresponding to transA2, cisA3, and transO1 from the transA1 conformation are presented. As suggested in a previous study,1 there appears to exist a significant potential barrier between the transA1 and cisA3 conformer, hindering conformational change between the two conformers. By contrast, the conformational change between the transA1 and transA2 appears to be easily accessible with temperature as the barrier height to the transition state from the transA1 toward the transA2 conformation is only 1.0 kcal/mol at the MP2/6-311++G** level. The effect of higher electron correlation beyond the B3LYP/6-31G* level on

the calculated barrier height appears to lower the barrier height as in the case of relative energies of the other conformers with respect to the energies of the transA1 conformer (see Table 2). For example, the barrier height to the transition state from the transA1 to the transA2 and cisA3 conformations changes from 1.9 and 28.6 kcal/mol at the B3LYP/6-31G* level to 1.0 and 25.3 kcal/mol, respectively, at the MP2/6-311++G ** level. These latter MP2 values of barrier heights were found to further decrease to 0.7 and 24.1 kcal/mol, respectively, at the CCSD(T)/ 6-311++G**//MP2/6-311++G** level calculations. In Table 5 we also presented the barrier heights to the transition states from the trans A1 to the transA2 and cisA3 conformations in aqueous solution phase computed under the PCM approximation19-24 at the B3LYP/6-31G* level. Interestingly, the barrier heights in the aqueous phase appear to be lower than the corresponding barrier heights in the gas phase in both cases, suggesting the possibility of rapid conformational changes between the transA1 and transA2 conformers in the aqueous solution phase at elevated temperatures. With regard to the population of the transO1 conformer compared to the population of more stable transA1 and transA2 conformers with temperature, one could expect from the results in Tables 4 and 5 that a nonnegligible portion of the transO1 conformer would be populated along with the most dominant transA1 conformer at elevated temperature as the Gibbs free energy difference between the transA1 and transO1 conformer at room temperature appears to be rather small (∼1 kcal/mol). To elucidate the rate of proton migration and possible proton migration pathway from the most stable transA1 form to the transO1 form of the protonated dialanine, we performed transition state optimizations connecting the transA1 and transO1 form along with the intrinsic reaction coordinate (IRC)19 calculation to confirm the validity of the obtained transition state at the B3LYP/6-31G* level. In Figure 3 the results of the IRC calculations are plotted along with the geometries of the reactant (transA1), product (transO1), and transition state. In the transition state the proton position appears to be close to the center between the O atom of the amide I CdO group and the N atom in the NH2 group. Similar transition state geometry was obtained at the MP2/6-311++G** level. The corresponding barrier height for proton migration from NH3+ to H+-O-C position appears to be rather high, 3.63 kcal/mol at the B3LYP/6-31G* level (2.9 kcal/mol at the MP2/6-311++G** level), suggesting a rather slow interconversion between the transA1 and transO1 conformations. Although this result is in general accord with the recent Car-Parrinello molecular dynamics study results which found rare interchange of proton between amino and carbonyl group in the MD trajectories of protonated dialanine near room temperature,2 high level ab initio calculations beyond DFT and MP2 levels appear necessary to clarify the mechanism of proton migration in this species. One of the fundamental issues related to structures and functions of peptides is concerned with the structural and energetic changes

7542

J. Phys. Chem. A, Vol. 114, No. 28, 2010

Sohn and Lee TABLE 6: Dihedral Angles (in degrees) and Relative Stability of the Protonated Dialanine Conformers in Watera transA1 Ψ1 ω Φo Ψ2 ∆Esol,relc ∆Eel,aqd ∆E298,aqe ∆G298,aqe

161.1 174.5 200.8 171.3 0 0 0 0

(-14.8)b (-2.8) (2.9) (-2.3)

transA2

transO1

cisA3

156.7 (-14.0) 167.5 (5.2) 282.8 (-5.9) 174.3 (-3.4) -0.73 -0.19 0.02(-0.17) -0.52(-0.71)

177.9 (-1.3) 178.5 (-0.7) 198.1 (2.4) 172.4 (-0.8) 5.18 6.42 -1.54(4.88) -1.05(5.37)

124.0 (-4.8) 336.9 (2.7) 302.1 (-3.8) 203.3 (-2.2) 1.81 1.94 0.31(2.25) 1.25(3.19)

a Based on the PCM results at the B3LYP/6-31G*level. b Values in parentheses represent the difference from the gas phase value. c ∆Esol,rel represents the relative solvation energy, which corresponds to the difference between the electronic energies in the gas and aqueous solution phase with respect to the salvation energy of the transA1 conformer. d ∆Eel,aq represents the relative electronic energy of the conformer with respect to the energy of the transA1 conformer in the aqueous solution (see the text). e ∆E298,aq and ∆G298,aq represent the relative thermal contribution to the energy and Gibbs free energy with respect to the transA1 conformer value at 298K and 1 atm. Values in parentheses represent the sum of thermal and electronic (∆Eel,aq) contribution.

by summing the ∆Eel values in Table 4 (reference relative electronic energies of conformers in the gas phase) and relative solvent effect ∆Esol,rel (difference between the electronic energies in th gas and aqueous solution phase with respect to the transA1 result) evaluated at the B3LYP/6-31G* level. From the ∆Eel,aq values in Table 6, the lowest energy conformer of the protonated dialanine in the aqueous environment appears to be the transA2 conformer. Both solvation (hydration) and thermal (entropic) effect appear to help stabilize the transA2 form, eventually making the transA2 form the lowest energy conformer in aqueous solution at room temperature. As a result, the transA2 conformer would be more populated than the transA1 conformer in aqueous solution at room temperature and the transA2 conformer would contribute to the spectral features of the protonated dialanine peptide ion more than the transA1 conformer in the aqueous solution phase, in contrast to the gas phase. Future experiment on the protonated dialanine in aqueous solution could check the premise. In the case of the transO1 and cisA3 conformers, hydration (transO1) and vibrational and thermal (cisA3) contributions appear to make these conformers lie in much higher energy states with respect to the energy of the transA1 (or transA2) conformer, practically making these states inaccessible at room temperature. This result appears to be in accord with the recent IRMPD experimental results on acetylcholine which found a reduction of the number of accessible conformers in aqueous solution compared to gas phase.39 Figure 3. IRC profile and optimized geometries of the transA1, transO1, and transition state (TS) at the B3LYP/6-31G* level.

IV. Summary and Conclusion

induced by the water solvent present in the aqueous solution state. To evaluate the effect of the water on the structure and stability of the protonated dialanine in aqueous solution, the PCM method19-24 was also employed. In Table 6 we present the changes in dihedral angles and relative stability between the conformers in aqueous solution, which was obtained from the calculations at the B3LYP/ 6-31G* level under the PCM approximation. In terms of structural changes, the results in Table 6 show that most dihedral angles do not change significantly except Ψ1 for transA1 and transA2 conformers, where a larger stabilization of energy is accompanied by solvation of the protonated dialanine in water compared to the transO1 or cisA3 conformer. To understand the energetic changes due to hydration correctly, one has to recognize that the relative electronic energies of conformers (∆Eel,aq) in Table 6 are obtained

We have examined the equilibrium structures and relative stability of the four lowest energy conformers (transA1, transA2, transO1, cisA3) of protonated dialanine cation in the gas phase using the HF, MP2, and hybrid DFT (BLYP, B3LYP, BPW91, B3PW91) methods as well as the potential barrier heights between transA1 and the other conformers and the effect of water solvent on the structures and relative stability of the four conformers employing the PCM method. The results of our work can be summarized as follows. (1) The conformer structure of the protonated dialanine, represented by four dihedral angles in this study, is sensitive to the electron correlation level and basis set employed for investigation. With respect to the reference dihedral angles which were taken as the MP2/aug-cc-pVTZ values, HF and hybrid DFT level results appear not to be close to the reference

Structure and Stability of Protonated Dialanine values within 3° and 2° (by average), respectively, regardless of the basis set employed. For the MP2 level optimization, use of a basis set at least with triple-ζ quality or inclusion of sufficient diffuse functions in the case of the double-ζ set appears necessary to obtain dihedral angles close to the reference values less than 2° by average. (2) In a similar vein, the accurate prediction of the relative stability of these conformers is strongly dependent upon the electron correlation degree and basis set quality. It was found after CCSD(T) correction that the lowest energetic transA1 and next to lowest, cisA3 conformer are almost isoenergetic, followed by transA2 (∼0.5 kcal/mol above transA1) and transO1 (∼1.2 kcal/mol above transA1) conformer in energy order. Considering computational feasibility along with reasonable accuracy, B3LYP/6-31G* appears to provide a good alternative to the more accurate (and expensive) method for the prediction of the structure and relative stability of protonated dialanine conformers. (3) The vibrational and thermal motions appear to be a critical factor in determining the relative stability of the conformers of the protonated dialanine at room temperature in the gas phase. After these nuclear motion contributions in the gas phase, transA1 becomes the lowest energy conformer, followed by transA2, transO1, and finally cisA3, in order. The small energy difference between transA1 and transA2 conformer suggests a significant population of transA2 conformer as well as transA1 conformer at room temperature, appearing to explain the nonsymmetric peak intensities due to CdO stretchings near 1750 cm-1 in the IRMPD spectrum of the protonated dialanine cation. (4) The potential barrier height between the most stable transA1 and transA2 conformers appears to be very small, probably less than 1 kcal/mol, suggesting protonated dialanine should exist as dynamical superposition of the transA1 and transA2 conformers at room temperature. By contrast, the barrier height between the transA1 and cisA3 conformers appears to be very high, probably more than 20 kcal/mol, indicating the interconversion between the two conformers would be rare and the population of the cisA3 conformer would be negligible near room temperature. This is in general accord with the apparent features of the IRMPD spectrum of the protonated dialanine. (5) The proton migration from the NH3+ group to the adjacent oxygen in the amide I carbonyl group would not be improbable at room temperature as the potential barrier height between the transA1 and transO1 conformer appears to be not significant, although the rate of migration would be slow as the barrier height appears to be much larger than thermal energy at room temperature. However, as this result is based on the B3LYP/6-31G* and MP2/ 6-311++G** calculations, high level ab initio calculations employing more accurate electron correlation treatment and larger basis set would be necessary to clarify this issue. (6) In aqueous solution, solvent effect (hydration), in combination with thermal effect, appears to significantly enhance the stability of the transA2 conformer, making the transA2 conformer more stable than the transA1 conformer at room temperature. As a result, transA2 becomes the most stable conformer, followed by transA1, cisA3, and finally transO1 in order. Structural change caused by hydration of protonated dialanine in water was not significant except for the transA1 and transA2 conformers for which changes of more than 10° in the dihedral angle defining NH3+ and CO group on the C-C skeleton were observed along with larger solvation energies compared to the transO1 or cisA3 conformer. As a result, the relative energies of the transO1 and cisA3 conformers with respect to the transA1 (or transA2) conformer become higher in aqueous solution compared to the gas phase.

J. Phys. Chem. A, Vol. 114, No. 28, 2010 7543 References and Notes (1) Lucas, B.; Gre´goire, G.; Lemaire, J.; Maıˆtre, P.; Ortega, J.-M.; Rupenyan, A.; Reimann, B.; Schermann, J. P.; Desfranc¸ois, C. Phys. Chem. Chem. Phys. 2004, 6, 2659. (2) Marinica, D. C.; Gre´goire, G.; Desfranc¸ois, C.; Schermann, J. P.; Borgis, D.; Gaigeot, M. P. J. Phys. Chem. A 2006, 110, 8802. (3) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (4) Stephens, J.; Devlin, F. J.; Chablowski, C. F.; Frisch, M. J. Phys. Chem. 1994, 98, 11623. (5) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (6) Becke, A. D. Phys. ReV. 1988, 38, 3098. (7) Moller, C.; Plesset, M. S. Phys. ReV. 1934, 46, 618. (8) Cizek, J. AdV. Chem. Phys. 1969, 14, 35. (9) Pople, J. A.; Head-Gordon, M.; Raghavachari, K. J. Chem. Phys. 1987, 87, 5968. (10) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. 1996, 117, 11225. (11) Maxwell, D. S.; Tirado-Rives, J.; Jorgensen, W. L. J. Comput. Chem. 1995, 16, 984. (12) Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. J. Comput. Chem 2005, 26, 1701. (13) Hess, B.; Kutzner, C.; Van Der Spoel, D.; Lindahl, E. J. Chem. Theory Comput. 2008, 4, 435. (14) Perdew, J. P.; Wang, Y. Phys. ReV. B 1992, 45, 13244. (15) Dunning, T. H., Jr. J. Chem. Phys. 1989, 90, 1007. (16) Woon, D. E.; Dunning, T. H., Jr. J. Chem. Phys. 1993, 98, 1358. (17) Pipano, A.; Gilman, R. R.; Shavitt, I. Chem. Phys. Lett. 1970, 5, 285. (18) Woon, D. E.; Dunning, T. H., Jr. J. Chem. Phys. 1995, 103, 4572. (19) Gonzalez, C.; Schelegel, H. B. J. Chem. Phys. 1989, 90, 2154. (20) Miertusˇ, S.; Scrocco, E.; Tomasi, J. Chem. Phys. 1981, 55, 117. (21) Pascual-Ahuir, J. L.; Silla, E.; Tun˜o´n, I. J. Comput. Chem. 1994, 15, 1127. (22) Mennucci, B.; Tomasi, J. J. Chem. Phys. 1997, 106, 5151. (23) Mennucci, B.; Cance`s, E.; Tomasi, J. J. Phys. Chem. B 1997, 101, 10506. (24) Miehlich, B.; Savin, A.; Stoll, H.; Preuss, H. Chem. Phys. Lett. 1989, 157, 200. (25) Cossi, M.; Barone, V.; Cammi, R.; Tomasi, J. Chem. Phys. Lett. 1996, 255, 327. (26) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, Revision D.01; Gaussian, Inc., Pittsburgh, PA, 2004. (27) Huh, S. B.; Lee, J. S. J. Chem. Phys. 2003, 118, 3035. (28) Truhlar, D. G. Chem. Phys. Lett. 1998, 98, 294. (29) Klopper, W.; Quack, M.; Suhm, M. A. Mol. Phys. 1998, 94, 105. (30) Klopper, W.; van Duijneveldt-van de Rigdt, J. G. C. M.; van Duijneveldt, F. B. Phys. Chem. Chem. Phys. 2000, 2, 2227. (31) Min, S. K.; Lee, E. C.; Lee, H. M.; Kim, D. Y.; Kim, D.; Kim, K. S. J. Comput. Chem. 2008, 29, 1208. (32) Johnson, R. D. In NIST Computational Chemistry Comparison and Benchmark Database; National Institute of Standards and Technology: Gaithersburg, MD, 2003. (33) Kohtani, M.; Schneider, J. E.; Jones, T. C.; Jarrold, M. F. J. Am. Chem. Soc. 2004, 126, 16981. (34) Wysocki, V. H.; Tsaprailis, G.; Smith, L. L.; Breci, L. A. J. Mass Spectrom. 2000, 35, 1399. (35) Paizs, B.; Csonka, I. P.; Lendvay, G.; Suhai, S. Rapid Commun. Mass Spectrom. 2001, 15, 637. (36) Schlosser, A.; Lehmann, W. D. J. Mass Spectrom. 2000, 35, 1382. (37) Johnson, R. S.; Martin, S. A.; Biemann, K. Int. J. Mass Spectrom. Ion Processes 1988, 86, 137. (38) Cox, K. A.; Gaskell, S. J.; Morris, M.; Whiting, A. J. Am. Soc. Mass Spectrom. 1996, 7, 522. (39) Seydou, M.; Gre´goire, G.; Liquier, J.; Lemaire, J.; Schermann, J. P.; Desfranc¸ois, C. J. Am. Chem. Soc. 2008, 130, 4187.

JP100380C