3476
J. Phys. Chem. B 2007, 111, 3476-3480
Are Isolated Nucleic Acid Bases Really Planar? A Car-Parrinello Molecular Dynamics Study Olexandr Isayev,† Al’ona Furmanchuk,† Oleg V. Shishkin,‡,⊥ Leond Gorb,†,§,⊥ and Jerzy Leszczynski*,†,⊥ Computational Center for Molecular Structure and Interactions, Jackson State UniVersity, Jackson, Mississippi 39217, STC “Institute for Single Crystals”, National Academy of Science of Ukraine, KharkiV 61001, Ukraine, Department of Molecular Biophysics, Institute of Molecular Biology and Genetics, National Academy of Science of Ukraine, KyiV 03143, Ukraine, and Ukrainian-American Laboratory of Computational Chemistry, KharkiV, Ukraine-Jackson, Mississippi 39217-0510 ReceiVed: January 31, 2007
Car-Parrinello molecular dynamics simulations of the flexibility of isolated DNA bases have been carried out. The comparison of lowest ring out-of-plane vibrations calculated by using MP2/cc-pvdz and BLYP/PW methods reveals that the DFT method with the plane wave basis set reasonably reproduces out-of-plane deformability of the pyrimidine ring in nucleic acid bases and could be used for reliable modeling of conformational flexibility of nucleobases. The conformational phase space of pyrimidine rings in thymine, cytosine, guanine, and adenine has been investigated by using the ab initio Car-Parrinello molecular dynamics method. It is demonstrated that all nucleic acid bases are highly flexible molecules and possess a nonplanar effective conformation of the pyrimidine ring despite the fact that the planar geometry corresponds to a minimum on the potential energy surface. The population of the planar geometry of the pyrimidine ring does not exceed 30%. Among the nonplanar conformations of the pyrimidine rings, the boat-like and half-chair conformations are the most populated.
Introduction The relationship between the structural dynamics of biomolecules and their biological function represents one of the key areas of research in modern structural biology.1 Knowledge of the rigid and flexible regions of a biomolecule allows detailed exploration of the ensemble of conformations available for a given state, and provides a means to predict changes in structural flexibility.2 It is also known that such fundamental biological processes like the replication of DNA and the transition between DNA-A and DNA-B are controlled by the conformational change of carbohydrate rings.1 Consequently, the flexibility of these biopolymers is associated with the relatively easy conformational transitions within the backbone such as rotation around σ bonds and transitions of sugar moiety between the south and north regions. Therefore, much effort has been spent relating the structure of DNA to sugar conformation.3 However, by using high-level ab initio calculations it was recently demonstrated that nucleic acid bases possess intrinsic conformational flexibility that results in easy out-of-plane deformation of the pyrimidine ring.4 It should be especially mentioned that, despite the fact that the described phenomenon is a typical dynamic effect, all referenced data have been obtained by analyzing the solutions of the time-independent Schro¨dinger equation. Therefore, they are not able to accurately predict the dynamical properties of nucleobases.5 * Address correspondence to this author. E-mail:
[email protected]. † Jackson State University. ‡ STC “Institute for Single Crystals”, National Academy of Science of Ukraine. § Department of Molecular Biophysics, Institute of Molecular Biology and Genetics, National Academy of Science of Ukraine. ⊥ Ukrainian-American Laboratory of Computational Chemistry.
From this perspective an accurate evaluation of nucleobase flexibility may be performed by molecular dynamic simulations. However, commonly used classic molecular dynamics methods strongly depend on parametrization of the force field.6 With the rapid progress of high performance computing the ab initio molecular dynamics simulations became accessible, in particular, the Car-Parrinello approach.7 This method has proven to be highly successful in biomodeling and materials science.8 Therefore our present goal is the complete characterization of the conformational phase space of nucleobases using ab initio molecular dynamics. We believe that the mapping of the key regions in the conformational phase space of nucleic acid bases is a necessary prerequisite for a deeper understanding of their biological functions. It also will be helpful for gleaning at least some hints about the possible motions of DNA molecules. Methods of Calculations The first principles molecular dynamics simulations within the Car-Parrinello scheme,7 as employed in the CPMD program,9 were carried out for the studied species. We used the BLYP functional, a plane wave (PW) basis set of TroullierMartins norm-conserving pseudopotentials;10 this approach will thereafter be referred to as BLYP/PW. The cluster boundary condition method of Martyna and Tuckerman11 was employed. The Kohn-Sham orbitals were expanded in a plane-wave basis set with a cutoff energy of 70 Ry. Both cutoff energy and the box sizes were adjusted for each system to yield energies converged to chemical accuracy. CP simulations were carried out at constant temperature with use of a massive Nose´-Hoover chains thermostat of length 4. It was checked that the fictitious electron kinetic energy remained small compared to that of the nuclei (adiabaticity) without the need for electronic thermostat-
10.1021/jp070857j CCC: $37.00 © 2007 American Chemical Society Published on Web 03/10/2007
Are Isolated Nucleic Acid Bases Really Planar?
J. Phys. Chem. B, Vol. 111, No. 13, 2007 3477
Figure 1. Representation of conformational space for the six-membered ring molecule, using the spherical coordinate {S, θ, Ψ} set. Basic ring conformations are indicated. The isolated segment of the sphere is the asymmetric unit.
ing (employing µ ) 700 amu). In all simulations, the time step was 4 au (about 0.097 fs). The total simulation time was about 10 ps for each species studied, including 2-3 ps of equilibration time. Conventional ab initio computations were carried out with the Gaussian 03 suite of programs.12 The geometrical parameters of all molecules were optimized by using the second-order Møller-Plesset perturbation theory13 with the standard cc-pvdz basis set14 basis sets. The vibrational frequencies were calculated within the harmonic approximation at the same level of theory. The conformation of the pyrimidine ring has been described in terms of puckering parameters developed by Zefirov, Palyulin, and Dashevskaya.15 There are three parameters S, θ, and Ψ characterizing six-membered-ring conformation (Figure 1). The S parameter describes the puckering degree of the ring. Its value is less than 0.1 if the values of endocyclic torsion angles do not exceed (5o. Therefore, conformations with S e 0.1 may be considered as planar or almost planar. Two polar angles θ and Ψ describe the type of ring conformation. In the case of some ideal conformations these angles adopt following values: chair, θ ) 0o, Ψ ) 0o; boat, θ ) 90o, Ψ ) 0o; twist-boat, θ ) 90o, Ψ ) 30o; sofa, θ ) 45o, Ψ ) 0o; and half-chair, θ ) 45o, Ψ ) 30o. However, it should be noted that these values of puckering parameters might be observed only in the case of ideal conformations. In real molecules one should consider some regions of values of polar angles around ideal values that should be assigned to specific conformation. Hence ranges of polar angle variations were taken as (15o for the θ angle and (10o for the Ψ angle. Conformations of ring with θ and Ψ values lying outside these regions were considered as intermediate conformations. However, more detailed analysis of such conformations indicates that they can be described as asymmetric analogues of half-chair, chair, and boat conformations. Therefore all these conformations may be summarized as, for example, boat-like or chair-like conformations. In the case of five-membered rings there are only two puckering parameters (S and Ψ) describing ring conformation. Like six-membered rings the S parameter describes the puckering degree of the ring. The Ψ polar angle reflects the variation of the five-membered ring conformation between two polar cases: envelope conformation (Ψ ) 0o) and twist conformation (Ψ ) 18o). Regions with ideal values of Ψ angle within (3o represent regions of ideal envelope and twist conformations. Other values of Ψ correspond to asymmetric twist conformation.
TABLE 1: Lowest Normal Modes, Zero Point Energy (ZPE), and Entropy of Nucleobases Calculated at MP2/ cc-pvdz and BLYP/PW Levels of Theory frequency (cm-1) MP2 PW thymine ν1 ν2 ν3 cytosine ν1 ν2 guanine ν1 ν2 ν3 adenine ν1 ν2
ZPE (au) assigment
MP2
PW
∆S (cal/(mol‚K)) MP2
PW
0.115605 0.112322 87.2
84.0
0.099404 0.096504 80.1
78.5
0.117611 0.113971 88.8
86.3
0.112882 0.109329 83.3
83.1
106 165 N(1) out-of-plane 136 253 N(3) out-of-plane 145 200 Me group rotation 129 240 N(1) out-of-plane 206 264 NH2 wag 136 202 CdO wag 157 232 N(1) out-of-plane 203 290 purine butterfly 163 159 NH2 wag 219 278 purine butterfly
Results and Discussion Validation of Computational Protocol with Respect to the Out-of-Plane Deformability of the Nucleobase. To verify the accuracy of the applied pseudopotential-plane waves approach utilized within the CPMD framework, it is necessary to compare the obtained results with the outcome of conventional, wellrecognized high-level quantum-chemical calculations. As such reference approximation we have chosen the second-order Møller-Plesset perturbation theory with the correlationconsistent cc-pvdz basis set. Therefore, the results obtained at the MP2/cc-pvdz level are compared with the similar BLYP/ PW ones. Keeping in mind the main goal of the present work we have compared the values of DNA bases’ lowest vibrational modes in the Table 1. These modes had been previously used to predict conformational flexibility and thermodynamical parameters of the DNA bases.4 We believe that they correctly describe the dynamic behavior of these molecular systems. An analysis of the data presented in Table 1 indicates reasonable agreement between BLYP/PW and MP2/cc-pvdz data. Both methods demonstrate that the pyrimidine ring outof-plane vibrations belong to the lowest vibrational mode of nucleobases. The largest differences are observed for thymine and cytosine. According to the BLYP/PW calculations the vibration corresponding to rotation of the methyl group has a
3478 J. Phys. Chem. B, Vol. 111, No. 13, 2007
Isayev et al.
Figure 2. Population of different values of the degree of puckering of the pyrimidine ring (a) and the imidazole ring (b) for corresponding nucleobases.
smaller frequency as compared to the N(3) out-of-plane vibration, contrary to MP2/cc-pvdz level predictions (Table 1). In the case of cytosine the two lowest vibrational modes calculated at the DFT and MP2 levels have an inverse order. It should be noted that all frequencies of out-of-plane and butterfly vibrations are systematically higher at the BLYP/PW level of theory as compared to the MP2 predictions and also to the known indirect experimental data (ref 16 and references therein). However, in electronic-structure calculations the quality of the results is also very sensitive to the particular basis set used. Converging the results obtained with different basis sets is not a trivial task for atomically localized orbitals. A plane-wave basis set has the advantage that the basis functions are naturally delocalized and orthogonal. The energy cutoff of the plane waves serves as a simple tuning parameter. Moreover, the absence of dependence of the plane waves on the atomic positions provides advantages in computing forces, vibrational frequencies, as well as intensities. Overall, application of the same methodology that we developed in ref 4 could result in prediction of more rigid DNA bases than was concluded by using the static ab initio approach. The ZPE and entropy values obtained within both approximations demonstrate reasonable agreement (Table 1). Differences in these values do not exceed 2.5 kcal/mol and 4 cal/(mol‚K) for ZPE and entropy values, respectively. Similar to vibrational frequencies these values are systematically smaller for BLYP/ PW calculations, indicating higher conformational rigidity of mucleobases at this level of theory. On the basis of this comparison, we can reasonably conclude that the BLYP/pseudopotential-plane waves approach is accurate enough to model the conformational flexibility of the pyrimidine ring. Therefore, the BLYP/PW method can be used for efficient energy evaluations during molecular dynamic simulation of nucleobases. Pyrimidine and Imidazole Ring Flexibility. CPMD simulation reveals that at room temperature all isolated nucleic acid bases are highly flexible molecules. The distribution of puckering parameter S for all DNA bases is depicted in Figure 2. The average degree of puckering of the pyrimidine ring is 0.140.20, and the corresponding values of the endocyclic torsion angles are up to (12-18o. Population of the planar geometry of this ring is only 12-30% (Table 2). Therefore, all nucleobases possess nonplanar effective conformations of the pyrimidine ring despite the planar geometry corresponding to a minimum on the potential energy surface.
TABLE 2: Population (%) of Different Conformations of Pyrimidine Rings in Nucleobases conformation
thymine
cytosine
guanine
adenine
planar symmetric boat asymmetric boat total boat-like total twist-boat symmetric chair asymmetric chair total chair-like symmetric half-chair asymmetric half-chair total half-chair-like sofa
11.9 6.0 6.0 12.0 29.0 2.5 5.0 7.5 10.8 17.2 28.0 11.6
17.1 5.7 5.6 11.3 24.2 3.6 3.9 7.5 13.5 16.0 29.5 10.4
17.8 5.5 5.6 11.1 24.0 4.0 4.7 8.7 10.7 17.3 28.0 10.4
30.3 5.7 5.2 10.9 24.4 2.4 3.2 5.6 9.1 11.1 20.2 8.6
An analysis of the nonplanar structures of the pyrimidine ring clearly demonstrates the presence of a very wide set of different ring conformations (Table 2). There are typical conformations of six-membered rings, like chair and boat, as well as numerous intermediate and distorted conformations, like sofa, twist-boat, half-chair, etc. Population analysis of individual pyrimidine ring conformers in nucleobases reveals that the twist-boat and asymmetric halfchair conformations are the most populated for all four molecules. It should be noted that these conformations together with sofa, symmetric boat, and half-chair conformations provide the main contribution to the total population of the nonplanar conformers of the pyrimidine ring. An analysis of the pathway of the described transformation indicates that these conformations are closely related (Figure 3a) because the transition between the two neighboring conformers requires minimal displacement of the atoms of the ring. This may be illustrated by using thymine as an example. Deviation of the N3 atom in the planar ring from the molecular plane results in the transition to a sofa conformation. Following this transition the atomic displacement occurs in the immediate proximity to the perturbed region of the ring thereby facilitating displacement of the second atom. There are two directions of displacement for the C2 atom. The first is characterized by the same direction as the N3 atom. This results in a twist-boat conformation. Otherwise, the opposite direction leads to an asymmetric half-chair conformation. When deviations of the C2 and N3 atoms from the mean plane of the remaining atoms of the ring become equal, the molecule adopts a symmetric boat or half-chair conformation. Therefore, one
Are Isolated Nucleic Acid Bases Really Planar?
J. Phys. Chem. B, Vol. 111, No. 13, 2007 3479
Figure 3. A scheme of conformational transitions of pyrimidine (a) and imidazole (b) rings from a planar equilibrium to different conformations.
TABLE 3: Population (%) of Different Conformations of Imidazole Rings in Guanine and Adenine conformation
guanine
adenine
planar symmetric twist asymmetric twist total twist-like envelope
66.0 9.1 14.9 24.0 10.0
64.2 8.8 16.7 25.5 10.3
can suggest that a high population of these conformers reflects the preferable direction of out-of-plane ring deformation. Taking into account the high conformational flexibility of the pyrimidine ring, it is interesting to analyze the behavior of the imidazole ring in guanine and adenine. Only three canonical conformations are possible for five-membered rings: planar and two nonplanar, namely twist and envelope (Figure 3b). Our results demonstrate that the imidazole ring is significantly more rigid that the pyrimidine ring. The range of variation of the degree of puckering is almost twice as small compared to that of the pyrimidine ring (Table 3). However, the average value of the degree of puckering (0.08) belongs to the region corresponding to the planar geometry with endocyclic torsion angles within a few degrees. Therefore one can conclude that the imidazole ring really possesses planar equilibrium geometry. The population values of the degree of puckering indicate that about 65% of the molecules possess a planar geometry (S < 0.1) (Table 3). However, a considerable number of structures with higher values of the degree of puckering are also observed that include various twist-like and envelope conformations. It should be noted that this distribution of populations is considerably more uniform compared to that of the pyrimidine ring. Conclusions We would like to point out that the ab initio molecular dynamics method coupled with the puckering analysis has made possible new insight into the conformational phase space of nucleabases and pathways of their conformation transformations. Our results clearly demonstrate that, despite the fact that the planar geometry corresponds to a minimum on the PES, the geometries of the isolated DNA bases are distributed between many conformations in some space under ambient conditions.
Again we would like to emphasize that this is purely a dynamical effect. Future work will examine the applicability and reliability of the current methodology and possible ways of extending it to more complex systems to account for the “real” DNA environment and bulk solvent effects. We expect that the conformational fluctuations lead to enhanced accessibility of the bases in the major and minor grooves of the helix compared with planar bases. In addition, high conformational flexibility of the pyrimidine ring represents one of the possible ways of geometry relaxation for the creation of the most favorable conditions for inter- and intramolecular interactions with other molecules like proteins, intercalators, drugs, etc. Acknowledgment. This work has been supported by NSF CREST (HRD-0318519), ONR (N00034-03-1-0116) grants, Mississippi Center for Supercomputer Research (MCRS), and AHPCRC (DAAH04-95-2-0003/contracts DAAH04-95-C-0008). Supporting Information Available: Short trajectory animation for one of the DNA bases (thymine) with projection onto the spherical coordinate {S, θ, Ψ} set. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) (a) Saenger, W. The Principles of Nucleic Acid Structure; Springer: Heidelberg, Germany, 1988. (b) Foloppe, N.; Nilsson, L.; MacKerell, A. F., Jr. Biopolymers 2002, 61, 61. (2) (a) Olson, W.; Gorin, A.; Lu, X.-J.; Hock, L.; Zhurkin, V. Proc. Natl. Acad. Sci. U.S.A. 1998, 95, 11163. (b) Strick, T. R.; Allemand, J.-F.; Bensimon, D.; Bensimon, A.; Croquette, V. Science 1996, 271, 1835. (3) (a) Prive´b, G.; Yanagi, K.; Dickerson, R. E. J. Mol. Biol. 1991, 217, 177. (b) Soliva, R.; Luque, F. J.; Alhambra, C.; Orozco, M. J. Biomol. Struct. Dyn. 1999, 17, 89. (4) (a) Shishkin, O. V. J. Chem. Soc., Chem. Commun. 1995, 1538. (b) Shishkin, O. V.; Gorb, L.; Hobza, P.; Leszczynski, J. Int. J. Quantum Chem. 2000, 80, 1116. (5) Shishkin, O.V.; Pichugin, K.; Gorb, L.; Leszczynski, J. J. Mol. Struct. 2002, 616 159. (6) (a) Cheatham, T. E., III; Young, M. A. Biopolymers 2000, 56, 232. (b) Feig, M.; Pettitt, B. M. Biophys. J. 1998, 75, 134. (7) Car, R.; Parrinello, M. Phys. ReV. Lett. 1985, 55, 2471. (8) (a) Carloni, P.; Rothlisberger, U.; Parrinello, M. Acc. Chem. Res. 2002, 35, 455. (b) Iftimie, R.; Minary, P.; Tuckerman, M. Proc. Natl. Acad. Sci. U.S.A., 2005, 102, 6654 and references cited therein.
3480 J. Phys. Chem. B, Vol. 111, No. 13, 2007 (9) CPMD v 3.9.1, Copyright IBM Corp 1990-2004; Copyright MPI fu¨r Festko¨rperforschung Stuttgart 1997-2001. (10) Troullier, N.; Martins, J. L. Phys. ReV. B 1991, 43, 1993. (11) Martyna, G. J.; Tuckerman, M. E. J. Chem. Phys. 1999, 110, 2810. (12) 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.; Bakken, V.; 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,
Isayev et al. 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 C.02; Gaussian, Inc.: Wallingford, CT, 2004. (13) Møller, C.; Plesset, M. S. Phys. ReV. 1934, 46, 618. (14) Woon, D. E.; Dunning, T. H. J. Chem. Phys. 1993, 98, 1358. (15) Zefirov, N. S.; Palyulin, V. A.; Dashevskaya, E. E. J. Phys. Org. Chem. 1990, 3, 147. (16) (a) Szczepaniak, K. S.; Szczesniak, M. M.; Person, W. B. J. Phys. Chem. A 2000, 104, 3852. (b) Aamouche, A.; Ghomi, M.; Colombeau, C.; Grajcar, L.; Baron, M. H.; Jobic, H.; Berthier, G. J. Phys. Chem. A 1997, 101, 1808. (c) Florian, J. J. Phys. Chem. 1993, 97, 10649.
