Stability and Interconversion of Acetylcholine Conformers

Abstract

The gas phase structures, energetics, and interconversion pathways of five lowest energy conformers of acetylcholine were examined employing the B3LYP, MP2, and CCSD(T) methods in conjunction with diverse basis sets including the correlation consistent aug-cc-pVDZ and aug-cc-pVTZ basis sets. It is found that use of adequate basis set containing proper polarization and diffuse functions capable of describing the floppy potential energy surface of acetylcholine is important in correctly predicting the relative stability of these conformers. The interconversion pathways and barrier heights between these conformers were elucidated by examining the potential energy surface for torsional motion, which also manifested the presence of chiral conformations of acetylcholine corresponding to the original conformations. On the basis of high level electronic energy calculations and thermal contribution analysis, four lowest energy conformers appear to be populated in the energy range of less than 1 kcal/mol at room temperature.

Introduction

Acetylcholine (2-Acetoxy-N,N,N-trimethylethanaminium) is an important neurotransmitter found in the peripheral and central nervous systems of vertebrate including humans1 and its conformational structure is crucial to understand the molecular recognition between the acetylcholine and its receptors in nervous system.2 It also serves as a mediator for parasympathetic function of the autonomic nerve system.3 Despite its importance and wide use in nervous system, however, the stable conformer structures of acetylcholine and their interconversion pathways, either in gas phase or in aqueous solution, were not well elucidated up to date.

This paper focuses on the accurate gas phase energetics and interconversion pathways of low energy conformers of positively charged acetylcholine ion using high level ab initio theoretical models. Although previous quantum theoretical studies4-11 indicated the presence of low energy zone depending upon orientation of acetoxy group in the system where stable conformations could be identified in the potential energy surface (PES) of acetylcholine, the energy differences between these conformations and interconversion mechanism could not be delineated as most of these previous studies employed the semiempirical and low level ab initio methods while the PES of acetylcholine is expected to be floppy. Another motivation for the present study is concerned with a recent experimental and theoretical study of the infrared multiphoton dissociation (IRMPD) spectrum of this ion, where analysis of spectral peaks was performed in terms of contributions from relevant low energy conformers.11

Previous theoretical studies on acetylcholine ion using DFT and MP2 methods9-11 have suggested that it could have four or five stable conformers within a few kilocalorie energy range and some of them appear to be close-lying, possibly within 1-2 kcal/mol. It has to be noted, however, that the energy difference of this small magnitude cannot be reliably predicted without proper account of electron correlation and basis set extension effect on relative energies between various conformations involved. Therefore, to determine the relative stability of low energy conformers accurately and to understand the conformational interconversion pathways between these close-lying conformers, we here employ more sophisticated theoretical models than used previously which must be able to manifest proper electron correlation and basis set extension effect on floppy potential energy surface of acetylcholine ion. In the present study we employ the MP212 and CCSD(T)13 methods as well as B3LYP14-17 method in conjunction with large basis sets to elucidate the energetics and interconversion pathways of close-lying conformers of this important species, which could also help us interpret the aforementioned IRMPD spectrum of acetylcholine correctly.

This paper is organized as follows. In the next section details of computational methods employed in this study are explained. In section III, the results of our computations and discussions about the results are presented, especially in relation to the spectral features of recently reported IRMPD spectrum of acetylcholine. Conclusion is in the final section.

 

Computational Details

To investigate the effect of electron correlation and basis set on the conformer structures and energetics of the aforementioned five lowest energy conformers (four of them considered to be within close energy range) of acetylcholine ion, three representative electron correlation methods, ab initio MP212 and CCSD(T) methods13 as well as B3LYP mrthod14-17 from density functional theory18 were employed in this study along with diverse basis sets including 6-31G*, 6-311++G** and correlation consistent basis sets aug-ccpVDZ and aug-cc-pVTZ.19,20 To examine the effect of diffuse functions on the optimized geometries of acetylcholine, we also employed the aug'-cc-pVDZ and aug'-cc-pVTZ sets which were constructed by removing the diffuse functions for hydrogen atom from the aug-cc-pVDZ and aug-cc-pVTZ sets, respectively. The largest basis set, aug-cc-pVTZ set, employed in this study is composed of [4s3p2d/5s4p3d2f] contracted functions.

After geometry optimization for each conformation, frequency calculations were performed to evaluate the thermal (Gibbs free energy) contributions to the electronic energy differences between various conformers. The frequency calculation results were also utilized to simulate the infrared spectrum of each conformer, which was compared with the experimental IRMPD spectrum for identification of contributions to the spectral structures of individual conformer. The transition states connecting various conformer structures were obtained by using various transition state optimization algorithms and verifying the first order saddle points of the transition structures through frequency calculations. To make sure the obtained transition structure is corresponding to the transition state connecting the desired initial and final conformer structures, intrinsic reaction coordinate (IRC) calculations21,22 were performed in both forward and reverse directions from the transition structure obtained. All correlated calculations were performed under the frozen core approximation. All calculations were performed using Gaussian program package.23

 

Results and Discussion

In Figure 1 the geometries of the five conformers of acetylcholine ion considered in this study (denoted C1 through C5) at the MP2/aug-cc-pVTZ optimization level are shown along with the dihedral angles θ1 (N(7)-C(6)-C(5)-O(4)) and θ2 (C(6)-C(5)-O(4)-C(2)), which represent the major structural features of acetylcholine as the acetoxy group is planar. In the case of θ1 for C2 and θ2 for C3 conformation, the MP2/aug-cc-pVTZ dihedral angles differ from the dihedral angles at the B3PW91/6-311++G** level11 by 10° and 7°, respectively, suggesting the order of stability of the conformers at the MP2/aug-cc-pVTZ level could be different from the results at the lower level of theory such as B3LYP/6-31G* or B3PW91/6-311++G**.11 This is well demonstrated in Table 1, where relative electronic energies (ΔE) of the other conformers with respect to C1 conformer at various levels of theory are presented.

Figure 1.Geometries of five conformers (C1-C5) of acetylcholine ion examined in this work. (θ1, θ2) are N(7)-C(6)-C(5)-O(4) and C(6)-C(5)-O(4)-C(2) torsion angles at the MP2/aug-cc-pVTZ level, respectively.

Table 1.aBasis sets contain diffuse functions for heavy atoms only except hydrogen. bMP2/aug-cc-pVTZ optimized geometries were used.

The first noticeable feature in Table 1 is that ΔE results are sensitive to the electron correlation method and basis set employed, especially to the quality of basis set. It is manifested that, regardless of electron correlation method adopted, employment of basis sets larger than the 6-311++G** set appears necessary to correctly predict the order of stability of these conformers of acetylcholine ion. It is worth noting that with large basis sets such as the aug-cc-pVDZ (or aug'-cc-pVDZ) and aug-cc-pVTZ (or aug'-cc-pVTZ), B3LYP, MP2, and CCSD(T) results all predict that the electronic energy of C5 is lower than that of C3, in contrast to the results with smaller basis sets. The energy difference between these two conformers (C3 and C5) appears to amount to 0.7 (MP2/aug-cc-pVTZ result)-1.0 (CCSD(T)/aug'-cc-pVDZ result) kcal/mol, although DFT results (B3LYP/aug-cc-pVDZ) exhibit smaller difference (~0.4 kcal/mol).

Table 2 shows more clearly how the relative stability of acetylcholine conformers can be affected by the size of basis set at the MP2 level. As the basis set becomes larger, the energy difference between C1 and C2 becomes larger, making C1 clearly the most stable conformation as well as reversing the order of stability between C3 and C5 conformers for this ion. It is interesting to note that the results with the aug'-cc-pVDZ set are similar to the results with the aug-cc-pVTZ set, suggesting the utility of the aug'-cc-pVDZ basis set in case use of more extensive basis set is not feasible.

Table 2.aGeometries optimized at the MP2/aug-cc-pVTZ are used throughout. bValues in parentheses represent the results at the geometries optimized at the MP2/6-311++G** level.

In Table 3 sum of Gibbs free energy contributions (computed at 298 K) to the reference electronic energies of each conformer, which were taken as the results at the CCSD(T)/aug'-cc-pVDZ level in Table 1, is presented. Although the thermal contribution results are different depending on the theoretical model employed, it is manifested that the energy of C3 conformer is higher than that of C1 (or C2 or C5), making C1, C2 and C5 most stable conformations for acetylcholine ion. The order of stability between these conformers (C1, C2 and C5) is not clear due to the strong dependence of thermal energy contribution results on the theoretical model employed, even though the large basis set calculations at the MP2 level appear to suggest C1 or C2 conformer as the lowest energy conformer for this ion.

Table 3.aElectronic energies were taken from the CCSD(T)/aug'-cc-pVDZ results in Table 1. bAll energies with respect to the energies for C1 conformer.

Given the energy differences between the various conformers of acetylcholine ion, it would be interesting to investigate the transition states and barrier heights connecting these conformations. For this purpose, transition state optimizations were performed at the MP2/6-311++G** level and the barrier heights were computed both at the MP2/6-311++G** and MP2/aug'-cc-pVDZ levels. The stationary points obtained were checked to yield only one imaginary frequency at the optimized geometries by frequency calculations. Finally, to make sure the optimized geometries obtained by this procedure are the desired transition state geometries connecting the two conformers of interest, intrinsic reaction coordinate calculations21,22 were carried out in forward and reverse directions beginning from the optimized transition state geometries.

In Table 4 and Figure 2 the barrier heights between various conformers containing C1, C2, C3, and C5 conformers and potential energy surface along conformational change pathways (θ1 and θ2) are presented. It is interesting to note from Figure 2 that there are no direct pathways connecting C1-C2 and C3-C5 transitions, for which C3 and C1 conformers play a role as an intermediate state, respectively. One important feature in Figure 2 is that, due to the symmetry of acetylcholine, there exist unsuperimposable mirror image conformers (represented as C1', C2', C3', and C5') for C1, C2, C3, and C5 conformers. The chirality of these two types of conformers is clearly seen in their IR and vibrational circular dichroism spectra in Figure 3. Thus, it is much easier for C2 conformer to change its conformation to C5' than C5. If the barrier for C2-C5' conformational change were significantly high for acetylcholine, these two types of conformers would be enantiomers separable at room temperature. However, the actual barrier height between C2-C5' conformational change in Table 4 suggests these two types of conformers would be difficult to separate at most temperatures of biological interest. Chirality in the absence of asymmetric carbon in the system has been known for other systems.24,25

Table 4.aBarrier heights with respect to the energies of the initial state conformers. bBarrier heights computed with the MP2/6-311++G** geometries. cValues in parentheses are the barrier heights with respect to the energies of the final state conformers.

Figure 2.Potential energy surface along θ1 and θ2 at the MP2/6-311++G** level. The energy points in this plot were obtained by optimizing all the geometric parameters except two dihedral angles (θ1, θ2). The contour interval is 0.5 kcal/mol.

Figure 3.(a) IR spectrum of C1 conformer (b) IR spectrum of C1' conformer (c) VCD spectrum of C1 conformer (d) VCD spectrum of C1' conformer. Y axis represents the peak intensity.

From the results in Table 3 and Table 4 it appears that C1, C2 and C5 conformers (compared to C3 and C4 conformers) would be dominantly populated at room temperature and IR spectrum of acetylcholine would exhibit the spectral features representing the contributions from these conformers. However, spectral contributions from C3 conformer also may not be negligible as the energy difference between C3 and C1 (or C2) appears to be small and the barrier heights from C3 to C1 or C2 appear to be higher than kT at room temperature. Therefore, it appears that the spectral features of IR spectrum for acetylcholine would have major contributions from these conformers.

Recently, Seydou et al. reported an IRMPD spectrum of acetylcholine ion between 800 cm−1 and 2000 cm−1 with assignment of peaks corresponding to stretching vibrations of backbone atoms in the acetylcholine conformers.11 The assignment was made by assuming C1, C2, and C3 conformers make most contributions to the spectral structures, which was based on molecular dynamics calculations at the HF/6-31G* level and simulation results of IR spectra of acetylcholine ion at the B3PW91/6-311++G** level. In view of the more accurate quantum chemical calculation results obtained in this work, however, it appears necessary to consider the contribution from C5 conformer as well as contributions from C1, C2, and C3 conformers to the spectral features in the IRMPD spectrum of acetylcholine ion for correct assignment of contributions from various conformers. For this purpose, we simulated the IR spectra of various conformers of acetylcholine ion at different levels of theory and identified the vibration frequencies having C(2)-O(3) and C(2)-O(4) stretching vibration as major contribution, respectively (see Fig. 1 for identification of relevant backbone atoms). Table 5 shows the aforementioned stretching frequencies at various levels of theory.

Table 5.aValues in parentheses represent the peak intensities.

In the case of C(2)-O(3) stretching vibration for which three Lorentzian shapes were used to fit the peaks near 1770 cm-1in the IRMPD spectrum (Fig. 4(c) in ref. 11), although the calculated absorption frequencies and intensities vary somewhat depending upon the level of theory employed, it appears that the conformers corresponding to the three peaks in the IRMPD spectrum in ref. 11 must be, from low to high frequency, C1(or C5), C2, and finally, C3 conformer, respectively, as the peak frequencies are shown to generally follow this order of calculated frequencies of aforementioned conformers. It is worth noting that C1 and C5 appear to have almost the same frequencies at the MP2 level (with larger basis sets), suggesting the lowest peak in the spectrum would correspond to contributions from both C1 and C5 conformers with larger peak intensity than the other remaining peaks. This is in good agreement with the spectral features in ref. 11 (Fig. 4(c)). Furthermore, the calculated intensity of C3 is stronger than the calculated intensity of C2, which also appears to be in good agreement with the spectral features in ref. 11, though the intensities in the IRMPD spectrum do not necessarily correspond to the absorption intensities in the IR spectrum due to the dependence of fragmentation yields on internal vibrational redistribution rates of the vibrational modes involved. Therefore, the spectral features around 1770 cm−1 in the IRMPD spectrum of acetylcholine ion appears to be in accord with our theoretical results about the stability of conformers of acetylcholine ion.

In the case of C(2)-O(4) stretching vibration, while the calculated frequencies in Table 5 appear to suggest more than two absorption peaks, the insufficient data points in the IRMPD spectrum (Fig. 4(b) in ref. 11) make it somewhat difficult to identify the accurate peak structures in this region. However, considering that the calculated frequencies of C1 and C5 conformer for C(2)-O(4) stretch are very close each other and the difference between the two peaks in Figure 4(b) in ref. 11 is similar to the difference in calculated frequencies of C1 (or C5) and C2 (or C3), it is suggested that the higher frequency peak should correspond to contributions from C1 and C5, with the lower frequency peak corresponding to contributions from C2 and C3. An IR spectrum of this ion with better resolution than the IRMPD spectrum in ref. 11 would be necessary to make a proper evaluation of the calculated results for C(2)-O(4) stretching vibration in Table 5.

 

Summary and Conclusion

We have examined the relative stability and interconversion of five lowest energy conformers (C1, C2, C3, C4, and C5) of acetylcholine ion in gas phase employing B3LYP, MP2, and CCSD(T) methods along with diverse basis sets including the correlation consistent aug-cc-pVDZ and augcc-pVTZ sets. Except for C4 lying in much higher energy region than the rest of the conformers, it was found that the energy order of the other conformers can be predicted correctly only when basis set is sufficiently large enough to contain a proper number of polarization and diffuse functions for correlated calculations. On the basis of MP2/aug-cc-pVTZ and CCSD(T)/aug'-cc-pVDZ results, the electronic energy of C5 conformer appears to be lower than the electronic energy of C3 by 0.5-1 kcal/mol, in contrast to the previous calculation results with smaller basis sets which suggested more stability for C3 than C5 conformer. The vibrational frequency and thermal contribution calculations also suggests that stability of C5 over C3 would be maintained at room temperature condition as the Gibbs free energy of C3 relative to C5 conformer at 298 K is shown to be larger by 0.3 to 0.8 kcal/mol based on at various levels of theory.

The presence of unsuperimposable configurations of mirror image for these conformers (C1', C2', C3', and C5') except C4 would make direct conversion from C1 to C2 (or from C3 to C5) difficult to proceed at room temperature without intermediate state of another conformer C3 (or C1). Similarly, while interconversion between C2 and C5 is much difficult to occur, C2 appears to change its conformation to C5' at room temperature without much difficulty due to the chirality of acetylcholine at room temperature.

The barrier heights between these conformers appear to be greater than kT at room temperature in most cases, suggesting distinct contributions to the IR spectral structure of acetylcholine ion from these conformers. This appears to be well manifested in the spectral features around 1770 cm−1 in the IRMPD spectrum of acetylcholine ion with which the general pattern of calculated frequencies and peak intensities of relevant conformers are in good accord. It has been suggested11 that presence of water molecules would reduce the number of accessible conformations in the case of acetylcholine. Since acetylcholine is known to be in hydrated form in neurotransmission of biological organism, it would be interesting to examine how the conformational energetics and interconversion pathways elucidated in this study would be affected by surrounding water molecules.