Atomistic understanding of the CT mismatched DNA base pair ...
WWW.C-CHEM.ORG
FULL PAPER
Atomistic Understanding of the CT Mismatched DNA Base Pair Tautomerization via the DPT: QM and QTAIM Computational Approaches Ol’ha O. Brovarets’,[a,b,c] and Dmytro M. Hovorun*[a,b,c] It was established that the cytosinethymine (CT) mismatched DNA base pair with cis-oriented N1H glycosidic bonds has propeller-like structure (jN3C4C4N3j 5 38.4 ), which is stabilized by three specific intermolecular interactions–two antiparallel N4H…O4 (5.19 kcal mol21) and N3H…N3 (6.33 kcal mol21) H-bonds and a van der Waals (vdW) contact O2…O2 (0.32 kcal mol21). The CT base mispair is thermodynamically stable structure (DGint 5 21.54 kcal mol21) and even slightly more stable than the AT Watson–Crick DNA base pair (DGint 5 21.43 kcal mol21) at the room temperature. It was shown that the CT $ C*T* tautomerization via the double proton transfer (DPT) is assisted by the O2…O2 vdW contact along the entire range of the intrinsic reaction coordinate (IRC). The positive value of the Grunenberg’s compliance constants ˚ /mdyn for the CT, C*T*, and (31.186, 30.265, and 22.166 A TSCT $ C*T*, respectively) proves that the O2…O2 vdW contact is a stabilizing interaction. Based on the sweeps of the H-bond energies, it was found that the N4H…O4/O4H…N4, and N3H…N3 H-bonds in the CT and C*T* base pairs are
anticooperative and weaken each other, whereas the middle N3H…N3 H-bond and the O2…O2 vdW contact are cooperative and mutually reinforce each other. It was found that the tautomerization of the CT base mispair through the DPT is concerted and asynchronous reaction that proceeds via the TSCT $ C*T* stabilized by the loosened N4AHAO4 covalent bridge, N3H…N3 H-bond (9.67 kcal mol21) and O2…O2 vdW contact (0.41 kcal mol21). The nine key points, describing the evolution of the CT $ C*T* tautomerization via the DPT, were detected and completely investigated along the IRC. The C*T* mispair was revealed to be the dynamically unstable structure with a lifetime 2.133 10213 s. In this case, as for the AT Watson–Crick DNA base pair, activates the mechanism of the quantum protection of the CT DNA base mispair from its spontaneous mutagenic tautomerization through the DPT. C 2013 Wiley Periodicals, Inc. V
Introduction
Without any exaggeration, a new era in structural studies of the nature of transitions and transversions began when the first X-ray data with acceptable resolution for short DNA fragments occured.[11,12] Since then, quite a number of structural and energetic data for all transitions and transversions in the double-stranded DNA without exception, which is in the crystalline state, as well as in the physiological aqueous solution, were obtained and analyzed by researchers using complementary physico-chemical methods.[13–23]
The nature of the formation of irregular pairs of nucleotide bases, that cause spontaneous point mutations in DNA,–the topic,[1–3] which continues to attract the attention of biologists, biochemists, and biophysicists nearly 60 years, starting since the time when Watson and Crick set the spatial architecture of DNA.[4] Freese was the first of researchers,[5] who proposed to classify all possible incorrect pairs of DNA bases as transitions and transversions. This classification has successfully stood the test of time and is now considered a classic. Current molecular genetic data suggest[6] that the total frequency of spontaneous point mutations in E.coli is 2.23 10210 per one DNA base pair that is synthesized per one generation. Thus, the frequency of transitions exceeds the frequency of transversions. At the same time, the frequency of spontaneous transversions is greater than the frequency of spontaneous transitions for the yeast Saccharomyces cerevisiae.[7] Interestingly, that among six possible transversions– three purinepurine [adenineadenine (AA), guanineguanine (GG), and adenineguanine (AG)] and three pyrimidinepyrimidine [cytosinecytosine (CC), thyminethymine (TT), and cytosinethymine (CT)]–GA, CT, and CC transversions are practically not deleted from the genome of E.coli.[8–10]
DOI: 10.1002/jcc.23412
[a] Ol’ha O. Brovarets’, D. M. Hovorun Department of Molecular and Quantum Biophysics, Institute of Molecular Biology and Genetics, National Academy of Sciences of Ukraine, 150 Zabolotnoho Str., Kyiv, 03680, Ukraine Email: [email protected] [b] Ol’ha O. Brovarets’, D. M. Hovorun Research and Educational Center, State Key Laboratory of Molecular and Cell Biology, 150 Zabolotnoho Str., Kyiv, 03680, Ukraine [c] Ol’ha O. Brovarets’, D. M. Hovorun Department of Molecular Biotechnology and Bioinformatics, Institute of High Technologies, Taras Shevchenko National University of Kyiv, 2 Hlushkova Ave., Kyiv, 03022, Ukraine Contract grant sponsor: State Fund for Fundamental Research of Ukraine (O.B.B.) (project no. GP/F56/074; Grant of the President of Ukraine); Contract grant sponsor: Ministry of Education and Science, Youth and Sports of Ukraine (Grant of the President of Ukraine) C 2013 Wiley Periodicals, Inc. V
Journal of Computational Chemistry 2013, DOI: 10.1002/jcc.23412
1
FULL PAPER
WWW.C-CHEM.ORG
Attempts were made to register experimentally all possible transitions and transversions in the base pair binding pocket of the DNA polymerase using X-ray analysis.[24,25] However, in contrast to the G*T and AC* purine–pyrimidine transitions,[26,27] for pyrimidine–pyrimidine transversions any structural data have been obtained in the base pair binding pocket of the DNA polymerase, which is in a closed conformation. The CT base pair, which are defined in the literature as short Watson–Crick base pair,[28] perhaps, arouses the discussions to a lesser extent among the pyrimidine–pyrimidine transversions.[14,16,17] Researchers converge on the fact that it is in the electroneutral state, is stabilized by two intermolecular N3H…N3 and N4H…O4 H-bonds (numbering of atoms is generally accepted[29]), and has the propeller-like structure in the composition of DNA under physiological conditions.[17] It is believed[14] that wedging of one water molecule into this pair from the side of the DNA minor groove and fixing it by two intermolecular H-bonds there adjusts its geometry to the Watson–Crick one, although direct experimental evidence is absent. It is clear that this result can not be transferred to a substantially hydrophobic base pair binding pocket of DNA polymerase, which is located in the replisome.[30–34] It is known that the CT pyrimidine–pyrimidine base mispair with cis-oriented N1H glycosidic bonds is the least stable among all three possible pyrimidine–pyrimidine base pairs in DNA.[14,16,17] In the recent work,[18] it was shown that anthracycline antibiotic nogalamycin stabilizes a mismatched CT base pair at the drug intercalation site in a DNA hairpin. At this, the CT base pair becomes planar. At the analysis of the possible point mutations that are generated by the CT transversion, it is traditionally believed that at the separation of the DNA strands and dissociation of this mispair into the monomers—C and T bases—their code properties remain unchanged. This approach completely ignores €wdin’s mechanism of the CT $ C*T* mutagenic tautothe Lo merization of this base mispair[35,36] (here and below the mutagenic tautomers[37–49]—imino in the case of C and enol in the case of T—are marked with asterisks). €wdin’s mechanism[35,36] Previously, we have shown that Lo for the AT Watson–Crick base pair did not pass the test of time.[45] This base pair does not tautomerize through the double proton transfer (DPT), because in this case a so-called quantum protection activates against the spontaneous mutagenic tautomerization.[44–46] However, as shows the experience in studying the mutagenic tautomerization mechanisms of the biologically important H-bonded pairs of nucleotide bases through the DPT,[46] these findings can not be automatically extended to the other base pair, even if they are close in their geometric and electronic structure. Thus, the triune purpose of our work, the results of which are set out below, is to study the structural and energetic properties of biologically important CT base mispair in the isolated state and the physico-chemical mechanism of its CT $ C*T* mutagenic tautomerization through the DPT. As the careful analysis of the literature shows, oddly enough, the results of this kind are missing, but they are very important, 2
Journal of Computational Chemistry 2013, DOI: 10.1002/jcc.23412
especially at the creating of new biosensors for identifying nucleobase mispairs in the double-stranded DNA.[48,50,51] Using the sweeps [the result of the scanning of the changes of the certain physico-chemical characteristic along the intrinsic reaction coordinate (IRC)] of the energetic, electrontopological, geometric and polar parameters of the CT $ C*T* tautomerization via the DPT along the IRC, we established that it is the concerted (i.e., this reaction involves no stable intermediates) and asynchronous (i.e., protons move with a time gap) process. If we take into account that replisoma most likely controls the canonical tautomeric status of €wdin’s mechanism of spontaneous the nucleobases,[49] the Lo point errors in DNA replication has no chances of being implemented in the CT base mispair.
Computational Methods All calculations have been carried out with the Gaussian’09 suite of programs.[52] Geometries and harmonic vibrational frequencies of the CT and C*T* short base pairs and the TSCT $ C*T* of their mutagenic tautomerization via the DPT were obtained using density functional theory (DFT)[53] with the B3LYP hybrid functional,[54] which includes Becke’s three-parameter exchange functional (B3)[55] combined with Lee, Yang, and Parr’s (LYP) correlation functional[56] in connection with Pople’s 6-31111G(d,p) basis set in vacuum. A scaling factor of 0.9668[39–44,46,57] has been used in the present work at the B3LYP quantummechanical (QM) level of theory to correct the harmonic frequencies of all the studied structures. We performed single point energy calculations at the correlated MP2 level of theory[58] with the 6-31111G(2df,pd) Pople’s[59–61] and cc-pVTZ/cc-pVQZ Dunning’s cc-type[62,63] basis sets for the B3LYP/6-31111G(d,p) geometries to consider electronic correlation effects as accurately as possible. MP2/6-31111G(2df,pd)//B3LYP/6-31111G(d,p), MP2/ cc-pVTZ//B3LYP/6-31111G(d,p), and MP2/cc-pVQZ//B3LYP/631111G(d,p) levels of theory have been successfully applied on similar systems recently studied and have been verified to give accurate normal mode frequencies, barrier heights, characteristics of intra- and intermolecular H-bonds, and geometries.[39–46,49,64–67] Moreover, an excellent agreement between computational and experimental NMR, UV, and IR spectroscopic data[37,57] evidences that the levels of theory applied for the single-point energy calculations (MP2/6-31111G(2df,pd), MP2/ cc-pVTZ, and MP2/cc-pVQZ), as well as the method used for the geometry optimization [B3LYP/6-31111G(d,p)] are reliable. The DFT method has been recommended in the literature for describing tautomerization phenomena of the H-bonded nucleobase pairs, as it has shown a good balance between computational cost and accuracy, and therefore, can be considered as the shortest way to MP2 results.[45,68,69] The correspondence of the stationary points to local minimum or TSCT $ C*T* on the potential energy landscape has been checked by the absence or the presence, respectively, of one and only one imaginary frequency corresponding to the normal mode that identifies the reaction coordinate. TSCT $ C*T* was located by means of synchronous transit-guided quasi-Newton method.[70,71] WWW.CHEMISTRYVIEWS.COM
FULL PAPER
WWW.C-CHEM.ORG
As the stationary points and TSCT $ C*T* were located, the reaction pathway was established by following the IRC in the forward and reverse directions from the TS using the Hessianbased predictor-corrector integration algorithm[72–74] with tight convergence criteria. These calculations eventually ensure that the proper reaction pathway, connecting the expected reactants and products on each side of the TS, has been found. We have investigated the evolution of the energetic, geometric, polar, and electron-topological characteristics of the H-bonds and base pairs along the reaction pathway establishing them at the each point of the IRC. The electronic interaction energies Eint have been computed at the MP2/6-31111G(2df,pd) level of theory for the geometries optimized at the DFT B3LYP/6-31111G(d,p) level of theory as the difference between the total energy of the base pair and the energies of the isolated monomers. In each case, the interaction energy was corrected for the basis set superposition error[75,76] through the counterpoise procedure.[77,78] The Gibbs free energy G values for all structures were obtained at room temperature (T 5 298.15 K) in the following way: G 5 Eel 1 Ecorr ;
(1)
where Eel-the electronic energy, whereas Ecorr-the thermal correction. The lifetime s of the tautomerized C*T* mispair can be estimated as 1/kr. The time s99.9% necessary to reach 99.9% of the equilibrium concentration of the CT reactant and the C*T* product of reaction in the system of reversible first-order forward (kf ) and reverse (kr) reactions was estimated by the formula[79]: s99:9% 5
ln103 : kf 1kr
(2)
To estimate the values of the forward kf and reverse kr rate constants for the CT$C*T* tautomerization reaction: kf;r 5C
kB T 2DDGf;r e RT h
(3)
we applied the standard TS theory,[79] in which quantum tunneling effects are accounted by the Wigner’s tunneling correction,[80] that is adequate for the DPT reactions[39–46,49]:
1 hmi C511 24 kB T;
program package AIMAll[82] with all the default options. Wave functions were obtained at the level of theory used for geometry optimization. The presence of a bond critical point (BCP),[81] namely the so-called (3,21) BCP and a bond path between hydrogen donor and acceptor, as well as the positive value of the Laplacian at this BCP (Dq 0), were considered as criteria for H-bond formation.[81] The energies of the conventional intermolecular H-bonds EHB in the CT and C*T* base pairs were evaluated by the empirical Iogansen’s formula[83]: EHB 50:33
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Dm240;
where Dm-the magnitude of the redshift (relative to the free molecule) of the stretching mode of the H-bonded groups involved in the H-bonding. The partial deuteration, namely the semideuteration of the amino group, was applied to eliminate the effect of vibrational resonances.[39–46,49,84] The energies of the O2…O2 van der Waals (vdW) contact EO2…O2 and of all intermolecular H-bonds EHB under the investigation of the sweeps of the their energies were evaluated by the empirical Espinosa–Molins–Lecomte (EML) formula[85,86] based on the electron density distribution at the (3,21) BCPs of the H-bonds and of the O2…O2 vdW contact: EHB=O2...O2 5 0:5 VðrÞ;
(4)
where kB-Boltzmann’s constant, T 5 298.15 K-temperature, h-Planck’s constant, DDGf,r-Gibbs free energy of activation for the forward and reverse DPT reactions (T 5 298.15 K), R-universal gas constant, mi-the magnitude of the imaginary frequency associated with the vibrational mode at the TS that connects reactants and products. Bader’s quantum theory “Atoms in molecules” (QTAIM) was applied to analyze electron density.[81] The topology of the electron density distribution has been examined in detail using
(6)
where V(r)-the value of a local potential energy density at the (3,21) BCPs. The energy of the N3H…N3 H-bond EN3H…N3 in the TSCT $ C*T* was estimated by the Nikolaienko–Bulavin–Hovorun formula[87]: EN3H...N3 5 22:03 1 225 q;
(7)
where q-the electron density at the (3,21) BCP of the N3H…N3 H-bond. Moreover, relative strength of the O2…O2 vdW contact was estimated by means of Grunenberg’s compliance constants formalism.[88,89] In contrast to force constants, the numerical values of compliance constants do not depend on the coordinate system. The physical meaning of compliance constants is deduced from their definition as partial second derivative of the potential energy due to an external force: Cij 5o2 E=ofi ofj :
2
(5)
(8)
In other words, compliance constants measure the displacement of an internal coordinate resulting from a unit force acting on it. As follows from this definition, a lower numerical value of compliance constant represents a stronger bond. The compliance constants were calculated using Compliance 3.0.2 program.[88,89] Period of the intermolecular vibrations T is calculated as 1 T5 mc , where m-the frequency of vibrations, c-the speed of the light in vacuum. The atomic numbering scheme for the C and T nucleobases is conventional.[29] Journal of Computational Chemistry 2013, DOI: 10.1002/jcc.23412
3
FULL PAPER
WWW.C-CHEM.ORG
Results and Discussion Structural and energetic peculiarities of the CT and C*T* base pairs and TS of their tautomerization via the DPT The discussion of the obtained results, which are presented in Tables (1-3) and Figures 1–9, we would like to start from the examination of the structural and energetic characteristics of the CT and C*T* base pairs and the transition state TSCT $ C*T* of the CT $ C*T* tautomerization. In this study, we consider only one of the two mirror-symmetric conformers (enantiomers) for each complex, as they are equivalent and do not require a separate analysis. The biologically important T and C*T* base mispairs, which are defined in the literature as short Watson–Crick base pairs,[28] are essentially nonplanar structures (/C4N3(T) N3C4(C) 5 38.4 for the T base pair and /C4N3(T*)N3C4(C*) 5 33.2 for the *T* base pair) with C1 symmetry and have propeller-like structures (Fig. 9). The heterocycles of the C and T bases as well as other DNA bases remain almost planar in the base pair, despite they are flexible molecules.[67,90,91] These base pairs are stabilized by the two antiparallel canonical Hbonds and a vdW contact: N4H…O4 (5.19 kcal mol21), N3H…N3 (6.33 kcal mol21), and O2…O2 (0.32 kcal mol21) for the T base pair and O4H…N4 (10.89 kcal mol21), N3H…N3 (6.34 kcal mol21), and O2…O2 (0.30 kcal mol21) for the *T* base pair (Tables 1 and 2). It is evident that the N3H…N3 Hbond is the strongest in the T base pair and the O4H…N4 H-bond is the strongest in the *T* base pair, whereas the O2…O2 vdW contact, exposed in the minor groove of the double-stranded DNA, is the weakest interaction in both complexes. The TSCT $ C*T* with imaginary frequency mi 5 725.9i cm21 has C1 symmetry (/C4N3(T)N3C4(C) 5 36.8 ) and is stabilized by the N4AHAO4 covalent bridge, conventional N3H…N3 H-bond (9.67 kcal mol21) and O2…O2 vdW contact (0.41 kcal mol21) (Tables 1 and 2, Fig. 1). The TSCT $ C*T* of the CT $ C*T* endothermic reaction is more similar to the C*T* base pair, obeying the Hammond’s postulate,[92] that states that endothermic reactions have a product-like TS.
The values of the electron density q at the BCPs of the Hbonds and vdW contacts, usually treated as a measure of these interactions,[84,93,94] fall within the range 0.002–0.066 a.u., that is consistent with the energetic data (Tables 1 and 2). The Laplacian of the electron density Dq is positive for all intrapair interactions and lies within the range 0.008–0.108 a.u. (Tables 1 and 2), demonstrating that H-bonds and vdW contact are attractive closed-shell interactions. The classical geometric criteria for the identifying of the Hbonds are satisfied for all AH…B H-bonds: the dH…B distances are less than the sum of Bondi’s vdW radii[95] of H (1.2 A˚) and ˚ ), N (1.55 A ˚ )] atoms. An elongation of the proton B [O (1.52 A donor group AH upon the formation of the conventional Hbond is positive in all cases and the angle of H-bonding is around 170 (for more details, see Tables 1 and 2). The spectroscopic data collected in Table 1 confirm geometric results. The shift in the frequency of the stretching mode of the AH donor group (the difference between the frequency for the AH group in monomer and in the base pair) is positive (shift to the red) for all H-bonds. It is interesting to note, that the T and C*T* base pairs examined in the present study are thermodynamically stable structures, as their Gibbs free energies of interaction (DGint 5 21.54 kcal mol21 (T) and DGint 5 29.16 kcal mol21 (*T*)) are less than zero, differing greatly from each other. Notably, the H-bonds and vdW contact make almost the same contribution to the stabilization of the T and C*T* base pairs [[REHB 1 EO2…O2]/jDEintj 5 85.4% (T) and 83.5% (*T*)]], whereas their electronic energies of the interaction differ drastically [DEint 5 213.86 kcal mol21 (T) and DEint 5 221.01 kcal mol21 (*T*)]. The energy relationship [REHB 1 EO2…O2]/jDEintj can not be considered as the unique physicochemical characteristic of the exceptionally T and C*T* base pairs and was also observed for the other H-bonded base pairs.[41,43–46,49,84] It was established by comparing the electronic energy of the back-barrier of the T $ C*T* tautomerization DDE 5 834.3 cm21 obtained at the MP2/cc-pVQZ//B3LYP/6-31111G(d,p) level of QM theory with the zero-point energy EZPE 5 1286.5 cm21 of
Table 1. Electron-topological, structural, vibrational, and energetic characteristics of the intermolecular H-bonds and O2. . .O2 vdW contact in the CT and C*T* base pairs and TSCT $ C*T* of their tautomerization obtained at the B3LYP/6-31111G(d,p) level of theory in vacuum.
Complex CT (IRC 5 29.32 Bohr)
TSCT
$ C*T*
(IRC 5 0.00 Bohr)
C*T* (IRC 5 3.70 Bohr)
AH…B -bond/ vdW contact
q[a]
Dq[b]
100e[c]
dA…B/dO2…O2[d]
dH…B[e]
DdAH[f ]
/AH…B[g]
Dm[h]
EHB/EO2…O2[i]
N4H…O4 N3H…N3 O2…O2 N3H…N3 O2…O2 O4H…N4 N3H…N3 O2…O2
0.031 0.028 0.002 0.052 0.002 0.066 0.029 0.002
0.108 0.078 0.008 0.092 0.009 0.096 0.081 0.008
3.48 6.50 21.46 5.95 17.67 4.69 6.32 26.46
2.873 3.001 3.730 2.781 3.639 2.637 2.974 3.744
1.851 1.974 – 1.728 – 1.611 1.959 –
0.016 0.024 – – – 0.060 0.024 –
174.8 170.1 – 167.6 – 173.9 165.9 –
287.3 408.2 – – – 1129.2 409.7 –
5.19 6.33 0.32* 9.67** 0.41* 10.89 6.34 0.30**
[a] The electron density at the BCP, a.u. [b] The Laplacian of the electron density at the BCP, a.u. [c] The ellipticity at the BCP. [d] Distance between A ˚ . [e] The distance between H and B (H-bond donor) and B (H-bond acceptor) atoms for H-bonds and between O2 oxygen atoms for the vdW contact, A atoms, A˚. [f ] The elongation of the H-bond donating group AH upon H-bonding, A˚. [g] The H-bond angle, degree. [h] The redshift of the stretching vibrational mode of the AH H-bonded and group, cm21. [i] Energy of the H-bonds, estimated by Iogansen’s formula or by Nikolaienko–Bulavin–Hovorun (marked with a double asterisk) formula energy of the O2…O2 vdW contact, estimated by Espinose–Molins–Lecomte formula (marked with an asterisk), kcal mol21.
4
Journal of Computational Chemistry 2013, DOI: 10.1002/jcc.23412
WWW.CHEMISTRYVIEWS.COM
FULL PAPER
WWW.C-CHEM.ORG
Table 2. Electron-topological, structural, and polar characteristics of the intermolecular H-bonds and O2. . .O2 vdW contact revealed in the structures of the nine key points obtained at the B3LYP/6-31111G(d,p) level of theory in vacuum.
Complex
AH…B -bond/ vdW contact
q
N4H…O4 N3…HN3 O2…O2 N4H…O4 N3…HN3 O2…O2 N4H…O4 N3-H-N3 O2…O2 N4H…O4 N3H…N3 O2…O2 N4H…O4 N3H…N3 O2…O2 N3H…N3 O2…O2 N4-H-O4 N3H…N3 O2…O2 N4…HO4 N3H…N3 O2…O2 N4…HO4 N3H…N3 O2…O2
0.031 0.028 0.002 0.065 0.108 0.003 0.067 0.144 0.003 0.071 0.103 0.003 0.134 0.053 0.002 0.052 0.002 0.163 0.051 0.002 0.117 0.048 0.002 0.066 0.029 0.002
Key point 1 (CT)(IRC 5 29.32 Bohr)
Key point 2 (DqN3…H 5 0, IRC 5 23.46 Bohr)
Key point 3 (qN3H 5 qHN3, IRC 5 23.24 Bohr)
Key point 4 (DqH…N3 5 0, IRC 5 22.98 Bohr)
Key point 5 (DqH…O4 5 0, IRC 5 20.06 Bohr)
Key point 6 (TSCT $
C*T*)
(IRC 5 0.00 Bohr)
Key point 7 (qN4H 5 qHO4, IRC 5 0.08 Bohr)
Key point 8 (DqN4…H 5 0, IRC 5 0.34 Bohr)
Key point 9 (C*T*) (IRC 5 3.70 Bohr)
Dq
100 E
dA…B/dO2…O2
dH…B
/AH…B
l
0.108 0.078 0.008 0.158 0.000 0.013 0.153 20.179 0.013 0.147 0.000 0.013 0.000 0.091 0.009 0.092 0.009 20.797 0.094 0.009 0.000 0.099 0.009 0.096 0.081 0.008
3.48 6.50 21.46 3.21 4.28 5.30 3.23 3.77 5.12 3.22 4.67 5.06 2.25 5.95 17.75 5.95 17.67 3.45 5.95 17.61 3.92 5.92 17.35 4.69 6.32 26.46
2.873 3.001 3.730 2.617 2.637 3.438 2.614 2.638 3.439 2.608 2.648 3.441 2.490 2.780 3.639 2.781 3.639 2.485 2.783 3.639 2.489 2.789 3.639 2.637 2.974 3.744
1.851 1.974 – 1.566 1.425 – 1.554 1.319 – 1.537 1.455 – 1.297 1.723 – 1.728 – 1.243 1.734 – 1.387 1.752 – 1.611 1.959 –
174.8 170.1 – 177.2 175.2 – 177.2 174.8 – 177.2 174.0 – 175.9 167.8 – 167.6 – 175.8 167.5 – 175.4 167.2 – 173.9 165.9 –
4.05
4.39
5.13
6.13
6.14
5.86 5.48
4.46
3.98
Notes: For footnote definitions see Table 1. l–the dipole moment of the complex, D.
the corresponding vibrational mode which frequency becomes imaginary in the TSCT $ C*T* that the C*T* mispair is dynamically (vibrationally) unstable[49,67,96] (Table 3, Fig. 2). Moreover, this statement is objective and does not depend on the chosed quantum-chemical level of theory (Table 3). It is obvious that these aforementioned facts put the possibility of the formation of the C*T* mispair involving mutagenic tautomers under a doubt. It should be also noted that the DDG 5 0.40 kcal mol21 Gibbs free energy of activation for the reverse reaction of the T $ C*T* tautomerization is less than kT 5 0.62 kcal mol21 in vacuum (DDG < kT). As a consequence (Table 3), the lifetime s of the C*T* mispair is 2.13 3 10213 s, that is much less than
the time, required for a replication machinery to forcibly dissociate a base pair into monomers (1029 ns[43]) during DNA replication. It means that the C*T* mispair “escapes form the hands” of the replication machinery due to its transformation to the CT base pair and further dissociation into the C and T monomers. Any of the six low-frequency intermolecular vibrations (23.9, 32.2, 57.5, 84.2, 99.1, and 102.8 cm21) can not develop during this lifetime, as the periods of these vibrations T (1.44 3 10212, 1.07 3 10212, 6.00 3 10213, 4.10 3 10213, 3.48 3 10213, 3.35 3 10213 s, respectively) are > the lifetime of the CT base pair s 5 2.13 3 10213 s. This observation additionally indicates that the C*T* base pair is dynamically unstable
Table 3. Energetic and kinetic characteristics of the CT $ C*T* tautomerization via the DPT in vacuo obtained at the different levels of QM theory. DDE[f ] Level of QM theory
DG[a] DE[b] DDGTSc DDETS[d] DDG[e] (kcal mol21) (cm21)
MP2/6-31111G(2df,pd)//B3LYP/6-31111G(d,p) 8.97 MP2/cc-pVTZ//B3LYP/6-31111G(d,p) 8.82 MP2/cc-pVQZ//B3LYP/6-31111G(d,p) 9.15
8.81 8.66 8.99
9.08 9.36 9.55
10.91 11.19 11.38
0.11 0.54 0.40
2.10 2.53 2.39
733.9 883.4 834.3
m[g]
EZPE[h]
s99.9%[j]
s[i] 213
2572.9 1286.5 1.3110 9.0610213 213 2572.9 1286.5 2.7010 1.8610212 2572.9 1286.5 2.1310213 1.4710212
[a] The relative Gibbs free energy of the C*T* base pair (DGCT 5 0; T 5 298.15 K), kcal mol21. [b] The relative electronic energy of the C*T* base pair (DECT 5 0), kcal mol21. [c] The Gibbs free energy of activation for the forward reaction of tautomerization (T 5 298.15 K), kcal mol21. [d] The activation electronic energy for the forward reaction of tautomerization, kcal mol21 [e] The Gibbs free energy of activation for the reverse reaction of tautomerization (T 5 298.15 K), kcal mol21. [f ] The activation electronic energy for the reverse reaction of tautomerization, kcal mol21. [g] The frequency of the vibrational mode in the C*T* base pair, which becomes imaginary in the TSCT $ C*T* of tautomerization obtained at the B3LYP/6-31111G(d,p) level of geometry optimization, cm21. [h] The zero-point vibrational energy associated with this normal mode, cm21. [i] The lifetime of the C*T* base pair, s. [j] The time necessary to reach 99.9% of the equilibrium concentration of the CT reactant and the C*T* product of the CT $ C*T* tautomerization reaction via the DPT, s.
Journal of Computational Chemistry 2013, DOI: 10.1002/jcc.23412
5
FULL PAPER
WWW.C-CHEM.ORG
Figure 1. Geometric structures of the nine key points describing the evolution of the CT $ C*T* tautomerization via the DPT along the IRC obtained at the B3LYP/6-31111G(d,p) level of theory in vacuo. Coordinates of each key point are presented below it. The dotted lines indicate AH…B H-bonds and O2…O2 vdW contact, whereas continuous lines show covalent bonds (their lengths are presented in angstroms). Carbon atoms are in light-blue, nitrogen in dark-blue, oxygen in red, and hydrogen in grey.
structure.[49,67,96] The time s99.9% necessary to reach 99.9% of the equilibrium concentration of the starting CT and the final C*T* base pair is equal to 1.47 3 10212 s (Table 3). The small value (
FULL PAPER
Atomistic Understanding of the CT Mismatched DNA Base Pair Tautomerization via the DPT: QM and QTAIM Computational Approaches Ol’ha O. Brovarets’,[a,b,c] and Dmytro M. Hovorun*[a,b,c] It was established that the cytosinethymine (CT) mismatched DNA base pair with cis-oriented N1H glycosidic bonds has propeller-like structure (jN3C4C4N3j 5 38.4 ), which is stabilized by three specific intermolecular interactions–two antiparallel N4H…O4 (5.19 kcal mol21) and N3H…N3 (6.33 kcal mol21) H-bonds and a van der Waals (vdW) contact O2…O2 (0.32 kcal mol21). The CT base mispair is thermodynamically stable structure (DGint 5 21.54 kcal mol21) and even slightly more stable than the AT Watson–Crick DNA base pair (DGint 5 21.43 kcal mol21) at the room temperature. It was shown that the CT $ C*T* tautomerization via the double proton transfer (DPT) is assisted by the O2…O2 vdW contact along the entire range of the intrinsic reaction coordinate (IRC). The positive value of the Grunenberg’s compliance constants ˚ /mdyn for the CT, C*T*, and (31.186, 30.265, and 22.166 A TSCT $ C*T*, respectively) proves that the O2…O2 vdW contact is a stabilizing interaction. Based on the sweeps of the H-bond energies, it was found that the N4H…O4/O4H…N4, and N3H…N3 H-bonds in the CT and C*T* base pairs are
anticooperative and weaken each other, whereas the middle N3H…N3 H-bond and the O2…O2 vdW contact are cooperative and mutually reinforce each other. It was found that the tautomerization of the CT base mispair through the DPT is concerted and asynchronous reaction that proceeds via the TSCT $ C*T* stabilized by the loosened N4AHAO4 covalent bridge, N3H…N3 H-bond (9.67 kcal mol21) and O2…O2 vdW contact (0.41 kcal mol21). The nine key points, describing the evolution of the CT $ C*T* tautomerization via the DPT, were detected and completely investigated along the IRC. The C*T* mispair was revealed to be the dynamically unstable structure with a lifetime 2.133 10213 s. In this case, as for the AT Watson–Crick DNA base pair, activates the mechanism of the quantum protection of the CT DNA base mispair from its spontaneous mutagenic tautomerization through the DPT. C 2013 Wiley Periodicals, Inc. V
Introduction
Without any exaggeration, a new era in structural studies of the nature of transitions and transversions began when the first X-ray data with acceptable resolution for short DNA fragments occured.[11,12] Since then, quite a number of structural and energetic data for all transitions and transversions in the double-stranded DNA without exception, which is in the crystalline state, as well as in the physiological aqueous solution, were obtained and analyzed by researchers using complementary physico-chemical methods.[13–23]
The nature of the formation of irregular pairs of nucleotide bases, that cause spontaneous point mutations in DNA,–the topic,[1–3] which continues to attract the attention of biologists, biochemists, and biophysicists nearly 60 years, starting since the time when Watson and Crick set the spatial architecture of DNA.[4] Freese was the first of researchers,[5] who proposed to classify all possible incorrect pairs of DNA bases as transitions and transversions. This classification has successfully stood the test of time and is now considered a classic. Current molecular genetic data suggest[6] that the total frequency of spontaneous point mutations in E.coli is 2.23 10210 per one DNA base pair that is synthesized per one generation. Thus, the frequency of transitions exceeds the frequency of transversions. At the same time, the frequency of spontaneous transversions is greater than the frequency of spontaneous transitions for the yeast Saccharomyces cerevisiae.[7] Interestingly, that among six possible transversions– three purinepurine [adenineadenine (AA), guanineguanine (GG), and adenineguanine (AG)] and three pyrimidinepyrimidine [cytosinecytosine (CC), thyminethymine (TT), and cytosinethymine (CT)]–GA, CT, and CC transversions are practically not deleted from the genome of E.coli.[8–10]
DOI: 10.1002/jcc.23412
[a] Ol’ha O. Brovarets’, D. M. Hovorun Department of Molecular and Quantum Biophysics, Institute of Molecular Biology and Genetics, National Academy of Sciences of Ukraine, 150 Zabolotnoho Str., Kyiv, 03680, Ukraine Email: [email protected] [b] Ol’ha O. Brovarets’, D. M. Hovorun Research and Educational Center, State Key Laboratory of Molecular and Cell Biology, 150 Zabolotnoho Str., Kyiv, 03680, Ukraine [c] Ol’ha O. Brovarets’, D. M. Hovorun Department of Molecular Biotechnology and Bioinformatics, Institute of High Technologies, Taras Shevchenko National University of Kyiv, 2 Hlushkova Ave., Kyiv, 03022, Ukraine Contract grant sponsor: State Fund for Fundamental Research of Ukraine (O.B.B.) (project no. GP/F56/074; Grant of the President of Ukraine); Contract grant sponsor: Ministry of Education and Science, Youth and Sports of Ukraine (Grant of the President of Ukraine) C 2013 Wiley Periodicals, Inc. V
Journal of Computational Chemistry 2013, DOI: 10.1002/jcc.23412
1
FULL PAPER
WWW.C-CHEM.ORG
Attempts were made to register experimentally all possible transitions and transversions in the base pair binding pocket of the DNA polymerase using X-ray analysis.[24,25] However, in contrast to the G*T and AC* purine–pyrimidine transitions,[26,27] for pyrimidine–pyrimidine transversions any structural data have been obtained in the base pair binding pocket of the DNA polymerase, which is in a closed conformation. The CT base pair, which are defined in the literature as short Watson–Crick base pair,[28] perhaps, arouses the discussions to a lesser extent among the pyrimidine–pyrimidine transversions.[14,16,17] Researchers converge on the fact that it is in the electroneutral state, is stabilized by two intermolecular N3H…N3 and N4H…O4 H-bonds (numbering of atoms is generally accepted[29]), and has the propeller-like structure in the composition of DNA under physiological conditions.[17] It is believed[14] that wedging of one water molecule into this pair from the side of the DNA minor groove and fixing it by two intermolecular H-bonds there adjusts its geometry to the Watson–Crick one, although direct experimental evidence is absent. It is clear that this result can not be transferred to a substantially hydrophobic base pair binding pocket of DNA polymerase, which is located in the replisome.[30–34] It is known that the CT pyrimidine–pyrimidine base mispair with cis-oriented N1H glycosidic bonds is the least stable among all three possible pyrimidine–pyrimidine base pairs in DNA.[14,16,17] In the recent work,[18] it was shown that anthracycline antibiotic nogalamycin stabilizes a mismatched CT base pair at the drug intercalation site in a DNA hairpin. At this, the CT base pair becomes planar. At the analysis of the possible point mutations that are generated by the CT transversion, it is traditionally believed that at the separation of the DNA strands and dissociation of this mispair into the monomers—C and T bases—their code properties remain unchanged. This approach completely ignores €wdin’s mechanism of the CT $ C*T* mutagenic tautothe Lo merization of this base mispair[35,36] (here and below the mutagenic tautomers[37–49]—imino in the case of C and enol in the case of T—are marked with asterisks). €wdin’s mechanism[35,36] Previously, we have shown that Lo for the AT Watson–Crick base pair did not pass the test of time.[45] This base pair does not tautomerize through the double proton transfer (DPT), because in this case a so-called quantum protection activates against the spontaneous mutagenic tautomerization.[44–46] However, as shows the experience in studying the mutagenic tautomerization mechanisms of the biologically important H-bonded pairs of nucleotide bases through the DPT,[46] these findings can not be automatically extended to the other base pair, even if they are close in their geometric and electronic structure. Thus, the triune purpose of our work, the results of which are set out below, is to study the structural and energetic properties of biologically important CT base mispair in the isolated state and the physico-chemical mechanism of its CT $ C*T* mutagenic tautomerization through the DPT. As the careful analysis of the literature shows, oddly enough, the results of this kind are missing, but they are very important, 2
Journal of Computational Chemistry 2013, DOI: 10.1002/jcc.23412
especially at the creating of new biosensors for identifying nucleobase mispairs in the double-stranded DNA.[48,50,51] Using the sweeps [the result of the scanning of the changes of the certain physico-chemical characteristic along the intrinsic reaction coordinate (IRC)] of the energetic, electrontopological, geometric and polar parameters of the CT $ C*T* tautomerization via the DPT along the IRC, we established that it is the concerted (i.e., this reaction involves no stable intermediates) and asynchronous (i.e., protons move with a time gap) process. If we take into account that replisoma most likely controls the canonical tautomeric status of €wdin’s mechanism of spontaneous the nucleobases,[49] the Lo point errors in DNA replication has no chances of being implemented in the CT base mispair.
Computational Methods All calculations have been carried out with the Gaussian’09 suite of programs.[52] Geometries and harmonic vibrational frequencies of the CT and C*T* short base pairs and the TSCT $ C*T* of their mutagenic tautomerization via the DPT were obtained using density functional theory (DFT)[53] with the B3LYP hybrid functional,[54] which includes Becke’s three-parameter exchange functional (B3)[55] combined with Lee, Yang, and Parr’s (LYP) correlation functional[56] in connection with Pople’s 6-31111G(d,p) basis set in vacuum. A scaling factor of 0.9668[39–44,46,57] has been used in the present work at the B3LYP quantummechanical (QM) level of theory to correct the harmonic frequencies of all the studied structures. We performed single point energy calculations at the correlated MP2 level of theory[58] with the 6-31111G(2df,pd) Pople’s[59–61] and cc-pVTZ/cc-pVQZ Dunning’s cc-type[62,63] basis sets for the B3LYP/6-31111G(d,p) geometries to consider electronic correlation effects as accurately as possible. MP2/6-31111G(2df,pd)//B3LYP/6-31111G(d,p), MP2/ cc-pVTZ//B3LYP/6-31111G(d,p), and MP2/cc-pVQZ//B3LYP/631111G(d,p) levels of theory have been successfully applied on similar systems recently studied and have been verified to give accurate normal mode frequencies, barrier heights, characteristics of intra- and intermolecular H-bonds, and geometries.[39–46,49,64–67] Moreover, an excellent agreement between computational and experimental NMR, UV, and IR spectroscopic data[37,57] evidences that the levels of theory applied for the single-point energy calculations (MP2/6-31111G(2df,pd), MP2/ cc-pVTZ, and MP2/cc-pVQZ), as well as the method used for the geometry optimization [B3LYP/6-31111G(d,p)] are reliable. The DFT method has been recommended in the literature for describing tautomerization phenomena of the H-bonded nucleobase pairs, as it has shown a good balance between computational cost and accuracy, and therefore, can be considered as the shortest way to MP2 results.[45,68,69] The correspondence of the stationary points to local minimum or TSCT $ C*T* on the potential energy landscape has been checked by the absence or the presence, respectively, of one and only one imaginary frequency corresponding to the normal mode that identifies the reaction coordinate. TSCT $ C*T* was located by means of synchronous transit-guided quasi-Newton method.[70,71] WWW.CHEMISTRYVIEWS.COM
FULL PAPER
WWW.C-CHEM.ORG
As the stationary points and TSCT $ C*T* were located, the reaction pathway was established by following the IRC in the forward and reverse directions from the TS using the Hessianbased predictor-corrector integration algorithm[72–74] with tight convergence criteria. These calculations eventually ensure that the proper reaction pathway, connecting the expected reactants and products on each side of the TS, has been found. We have investigated the evolution of the energetic, geometric, polar, and electron-topological characteristics of the H-bonds and base pairs along the reaction pathway establishing them at the each point of the IRC. The electronic interaction energies Eint have been computed at the MP2/6-31111G(2df,pd) level of theory for the geometries optimized at the DFT B3LYP/6-31111G(d,p) level of theory as the difference between the total energy of the base pair and the energies of the isolated monomers. In each case, the interaction energy was corrected for the basis set superposition error[75,76] through the counterpoise procedure.[77,78] The Gibbs free energy G values for all structures were obtained at room temperature (T 5 298.15 K) in the following way: G 5 Eel 1 Ecorr ;
(1)
where Eel-the electronic energy, whereas Ecorr-the thermal correction. The lifetime s of the tautomerized C*T* mispair can be estimated as 1/kr. The time s99.9% necessary to reach 99.9% of the equilibrium concentration of the CT reactant and the C*T* product of reaction in the system of reversible first-order forward (kf ) and reverse (kr) reactions was estimated by the formula[79]: s99:9% 5
ln103 : kf 1kr
(2)
To estimate the values of the forward kf and reverse kr rate constants for the CT$C*T* tautomerization reaction: kf;r 5C
kB T 2DDGf;r e RT h
(3)
we applied the standard TS theory,[79] in which quantum tunneling effects are accounted by the Wigner’s tunneling correction,[80] that is adequate for the DPT reactions[39–46,49]:
1 hmi C511 24 kB T;
program package AIMAll[82] with all the default options. Wave functions were obtained at the level of theory used for geometry optimization. The presence of a bond critical point (BCP),[81] namely the so-called (3,21) BCP and a bond path between hydrogen donor and acceptor, as well as the positive value of the Laplacian at this BCP (Dq 0), were considered as criteria for H-bond formation.[81] The energies of the conventional intermolecular H-bonds EHB in the CT and C*T* base pairs were evaluated by the empirical Iogansen’s formula[83]: EHB 50:33
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Dm240;
where Dm-the magnitude of the redshift (relative to the free molecule) of the stretching mode of the H-bonded groups involved in the H-bonding. The partial deuteration, namely the semideuteration of the amino group, was applied to eliminate the effect of vibrational resonances.[39–46,49,84] The energies of the O2…O2 van der Waals (vdW) contact EO2…O2 and of all intermolecular H-bonds EHB under the investigation of the sweeps of the their energies were evaluated by the empirical Espinosa–Molins–Lecomte (EML) formula[85,86] based on the electron density distribution at the (3,21) BCPs of the H-bonds and of the O2…O2 vdW contact: EHB=O2...O2 5 0:5 VðrÞ;
(4)
where kB-Boltzmann’s constant, T 5 298.15 K-temperature, h-Planck’s constant, DDGf,r-Gibbs free energy of activation for the forward and reverse DPT reactions (T 5 298.15 K), R-universal gas constant, mi-the magnitude of the imaginary frequency associated with the vibrational mode at the TS that connects reactants and products. Bader’s quantum theory “Atoms in molecules” (QTAIM) was applied to analyze electron density.[81] The topology of the electron density distribution has been examined in detail using
(6)
where V(r)-the value of a local potential energy density at the (3,21) BCPs. The energy of the N3H…N3 H-bond EN3H…N3 in the TSCT $ C*T* was estimated by the Nikolaienko–Bulavin–Hovorun formula[87]: EN3H...N3 5 22:03 1 225 q;
(7)
where q-the electron density at the (3,21) BCP of the N3H…N3 H-bond. Moreover, relative strength of the O2…O2 vdW contact was estimated by means of Grunenberg’s compliance constants formalism.[88,89] In contrast to force constants, the numerical values of compliance constants do not depend on the coordinate system. The physical meaning of compliance constants is deduced from their definition as partial second derivative of the potential energy due to an external force: Cij 5o2 E=ofi ofj :
2
(5)
(8)
In other words, compliance constants measure the displacement of an internal coordinate resulting from a unit force acting on it. As follows from this definition, a lower numerical value of compliance constant represents a stronger bond. The compliance constants were calculated using Compliance 3.0.2 program.[88,89] Period of the intermolecular vibrations T is calculated as 1 T5 mc , where m-the frequency of vibrations, c-the speed of the light in vacuum. The atomic numbering scheme for the C and T nucleobases is conventional.[29] Journal of Computational Chemistry 2013, DOI: 10.1002/jcc.23412
3
FULL PAPER
WWW.C-CHEM.ORG
Results and Discussion Structural and energetic peculiarities of the CT and C*T* base pairs and TS of their tautomerization via the DPT The discussion of the obtained results, which are presented in Tables (1-3) and Figures 1–9, we would like to start from the examination of the structural and energetic characteristics of the CT and C*T* base pairs and the transition state TSCT $ C*T* of the CT $ C*T* tautomerization. In this study, we consider only one of the two mirror-symmetric conformers (enantiomers) for each complex, as they are equivalent and do not require a separate analysis. The biologically important T and C*T* base mispairs, which are defined in the literature as short Watson–Crick base pairs,[28] are essentially nonplanar structures (/C4N3(T) N3C4(C) 5 38.4 for the T base pair and /C4N3(T*)N3C4(C*) 5 33.2 for the *T* base pair) with C1 symmetry and have propeller-like structures (Fig. 9). The heterocycles of the C and T bases as well as other DNA bases remain almost planar in the base pair, despite they are flexible molecules.[67,90,91] These base pairs are stabilized by the two antiparallel canonical Hbonds and a vdW contact: N4H…O4 (5.19 kcal mol21), N3H…N3 (6.33 kcal mol21), and O2…O2 (0.32 kcal mol21) for the T base pair and O4H…N4 (10.89 kcal mol21), N3H…N3 (6.34 kcal mol21), and O2…O2 (0.30 kcal mol21) for the *T* base pair (Tables 1 and 2). It is evident that the N3H…N3 Hbond is the strongest in the T base pair and the O4H…N4 H-bond is the strongest in the *T* base pair, whereas the O2…O2 vdW contact, exposed in the minor groove of the double-stranded DNA, is the weakest interaction in both complexes. The TSCT $ C*T* with imaginary frequency mi 5 725.9i cm21 has C1 symmetry (/C4N3(T)N3C4(C) 5 36.8 ) and is stabilized by the N4AHAO4 covalent bridge, conventional N3H…N3 H-bond (9.67 kcal mol21) and O2…O2 vdW contact (0.41 kcal mol21) (Tables 1 and 2, Fig. 1). The TSCT $ C*T* of the CT $ C*T* endothermic reaction is more similar to the C*T* base pair, obeying the Hammond’s postulate,[92] that states that endothermic reactions have a product-like TS.
The values of the electron density q at the BCPs of the Hbonds and vdW contacts, usually treated as a measure of these interactions,[84,93,94] fall within the range 0.002–0.066 a.u., that is consistent with the energetic data (Tables 1 and 2). The Laplacian of the electron density Dq is positive for all intrapair interactions and lies within the range 0.008–0.108 a.u. (Tables 1 and 2), demonstrating that H-bonds and vdW contact are attractive closed-shell interactions. The classical geometric criteria for the identifying of the Hbonds are satisfied for all AH…B H-bonds: the dH…B distances are less than the sum of Bondi’s vdW radii[95] of H (1.2 A˚) and ˚ ), N (1.55 A ˚ )] atoms. An elongation of the proton B [O (1.52 A donor group AH upon the formation of the conventional Hbond is positive in all cases and the angle of H-bonding is around 170 (for more details, see Tables 1 and 2). The spectroscopic data collected in Table 1 confirm geometric results. The shift in the frequency of the stretching mode of the AH donor group (the difference between the frequency for the AH group in monomer and in the base pair) is positive (shift to the red) for all H-bonds. It is interesting to note, that the T and C*T* base pairs examined in the present study are thermodynamically stable structures, as their Gibbs free energies of interaction (DGint 5 21.54 kcal mol21 (T) and DGint 5 29.16 kcal mol21 (*T*)) are less than zero, differing greatly from each other. Notably, the H-bonds and vdW contact make almost the same contribution to the stabilization of the T and C*T* base pairs [[REHB 1 EO2…O2]/jDEintj 5 85.4% (T) and 83.5% (*T*)]], whereas their electronic energies of the interaction differ drastically [DEint 5 213.86 kcal mol21 (T) and DEint 5 221.01 kcal mol21 (*T*)]. The energy relationship [REHB 1 EO2…O2]/jDEintj can not be considered as the unique physicochemical characteristic of the exceptionally T and C*T* base pairs and was also observed for the other H-bonded base pairs.[41,43–46,49,84] It was established by comparing the electronic energy of the back-barrier of the T $ C*T* tautomerization DDE 5 834.3 cm21 obtained at the MP2/cc-pVQZ//B3LYP/6-31111G(d,p) level of QM theory with the zero-point energy EZPE 5 1286.5 cm21 of
Table 1. Electron-topological, structural, vibrational, and energetic characteristics of the intermolecular H-bonds and O2. . .O2 vdW contact in the CT and C*T* base pairs and TSCT $ C*T* of their tautomerization obtained at the B3LYP/6-31111G(d,p) level of theory in vacuum.
Complex CT (IRC 5 29.32 Bohr)
TSCT
$ C*T*
(IRC 5 0.00 Bohr)
C*T* (IRC 5 3.70 Bohr)
AH…B -bond/ vdW contact
q[a]
Dq[b]
100e[c]
dA…B/dO2…O2[d]
dH…B[e]
DdAH[f ]
/AH…B[g]
Dm[h]
EHB/EO2…O2[i]
N4H…O4 N3H…N3 O2…O2 N3H…N3 O2…O2 O4H…N4 N3H…N3 O2…O2
0.031 0.028 0.002 0.052 0.002 0.066 0.029 0.002
0.108 0.078 0.008 0.092 0.009 0.096 0.081 0.008
3.48 6.50 21.46 5.95 17.67 4.69 6.32 26.46
2.873 3.001 3.730 2.781 3.639 2.637 2.974 3.744
1.851 1.974 – 1.728 – 1.611 1.959 –
0.016 0.024 – – – 0.060 0.024 –
174.8 170.1 – 167.6 – 173.9 165.9 –
287.3 408.2 – – – 1129.2 409.7 –
5.19 6.33 0.32* 9.67** 0.41* 10.89 6.34 0.30**
[a] The electron density at the BCP, a.u. [b] The Laplacian of the electron density at the BCP, a.u. [c] The ellipticity at the BCP. [d] Distance between A ˚ . [e] The distance between H and B (H-bond donor) and B (H-bond acceptor) atoms for H-bonds and between O2 oxygen atoms for the vdW contact, A atoms, A˚. [f ] The elongation of the H-bond donating group AH upon H-bonding, A˚. [g] The H-bond angle, degree. [h] The redshift of the stretching vibrational mode of the AH H-bonded and group, cm21. [i] Energy of the H-bonds, estimated by Iogansen’s formula or by Nikolaienko–Bulavin–Hovorun (marked with a double asterisk) formula energy of the O2…O2 vdW contact, estimated by Espinose–Molins–Lecomte formula (marked with an asterisk), kcal mol21.
4
Journal of Computational Chemistry 2013, DOI: 10.1002/jcc.23412
WWW.CHEMISTRYVIEWS.COM
FULL PAPER
WWW.C-CHEM.ORG
Table 2. Electron-topological, structural, and polar characteristics of the intermolecular H-bonds and O2. . .O2 vdW contact revealed in the structures of the nine key points obtained at the B3LYP/6-31111G(d,p) level of theory in vacuum.
Complex
AH…B -bond/ vdW contact
q
N4H…O4 N3…HN3 O2…O2 N4H…O4 N3…HN3 O2…O2 N4H…O4 N3-H-N3 O2…O2 N4H…O4 N3H…N3 O2…O2 N4H…O4 N3H…N3 O2…O2 N3H…N3 O2…O2 N4-H-O4 N3H…N3 O2…O2 N4…HO4 N3H…N3 O2…O2 N4…HO4 N3H…N3 O2…O2
0.031 0.028 0.002 0.065 0.108 0.003 0.067 0.144 0.003 0.071 0.103 0.003 0.134 0.053 0.002 0.052 0.002 0.163 0.051 0.002 0.117 0.048 0.002 0.066 0.029 0.002
Key point 1 (CT)(IRC 5 29.32 Bohr)
Key point 2 (DqN3…H 5 0, IRC 5 23.46 Bohr)
Key point 3 (qN3H 5 qHN3, IRC 5 23.24 Bohr)
Key point 4 (DqH…N3 5 0, IRC 5 22.98 Bohr)
Key point 5 (DqH…O4 5 0, IRC 5 20.06 Bohr)
Key point 6 (TSCT $
C*T*)
(IRC 5 0.00 Bohr)
Key point 7 (qN4H 5 qHO4, IRC 5 0.08 Bohr)
Key point 8 (DqN4…H 5 0, IRC 5 0.34 Bohr)
Key point 9 (C*T*) (IRC 5 3.70 Bohr)
Dq
100 E
dA…B/dO2…O2
dH…B
/AH…B
l
0.108 0.078 0.008 0.158 0.000 0.013 0.153 20.179 0.013 0.147 0.000 0.013 0.000 0.091 0.009 0.092 0.009 20.797 0.094 0.009 0.000 0.099 0.009 0.096 0.081 0.008
3.48 6.50 21.46 3.21 4.28 5.30 3.23 3.77 5.12 3.22 4.67 5.06 2.25 5.95 17.75 5.95 17.67 3.45 5.95 17.61 3.92 5.92 17.35 4.69 6.32 26.46
2.873 3.001 3.730 2.617 2.637 3.438 2.614 2.638 3.439 2.608 2.648 3.441 2.490 2.780 3.639 2.781 3.639 2.485 2.783 3.639 2.489 2.789 3.639 2.637 2.974 3.744
1.851 1.974 – 1.566 1.425 – 1.554 1.319 – 1.537 1.455 – 1.297 1.723 – 1.728 – 1.243 1.734 – 1.387 1.752 – 1.611 1.959 –
174.8 170.1 – 177.2 175.2 – 177.2 174.8 – 177.2 174.0 – 175.9 167.8 – 167.6 – 175.8 167.5 – 175.4 167.2 – 173.9 165.9 –
4.05
4.39
5.13
6.13
6.14
5.86 5.48
4.46
3.98
Notes: For footnote definitions see Table 1. l–the dipole moment of the complex, D.
the corresponding vibrational mode which frequency becomes imaginary in the TSCT $ C*T* that the C*T* mispair is dynamically (vibrationally) unstable[49,67,96] (Table 3, Fig. 2). Moreover, this statement is objective and does not depend on the chosed quantum-chemical level of theory (Table 3). It is obvious that these aforementioned facts put the possibility of the formation of the C*T* mispair involving mutagenic tautomers under a doubt. It should be also noted that the DDG 5 0.40 kcal mol21 Gibbs free energy of activation for the reverse reaction of the T $ C*T* tautomerization is less than kT 5 0.62 kcal mol21 in vacuum (DDG < kT). As a consequence (Table 3), the lifetime s of the C*T* mispair is 2.13 3 10213 s, that is much less than
the time, required for a replication machinery to forcibly dissociate a base pair into monomers (1029 ns[43]) during DNA replication. It means that the C*T* mispair “escapes form the hands” of the replication machinery due to its transformation to the CT base pair and further dissociation into the C and T monomers. Any of the six low-frequency intermolecular vibrations (23.9, 32.2, 57.5, 84.2, 99.1, and 102.8 cm21) can not develop during this lifetime, as the periods of these vibrations T (1.44 3 10212, 1.07 3 10212, 6.00 3 10213, 4.10 3 10213, 3.48 3 10213, 3.35 3 10213 s, respectively) are > the lifetime of the CT base pair s 5 2.13 3 10213 s. This observation additionally indicates that the C*T* base pair is dynamically unstable
Table 3. Energetic and kinetic characteristics of the CT $ C*T* tautomerization via the DPT in vacuo obtained at the different levels of QM theory. DDE[f ] Level of QM theory
DG[a] DE[b] DDGTSc DDETS[d] DDG[e] (kcal mol21) (cm21)
MP2/6-31111G(2df,pd)//B3LYP/6-31111G(d,p) 8.97 MP2/cc-pVTZ//B3LYP/6-31111G(d,p) 8.82 MP2/cc-pVQZ//B3LYP/6-31111G(d,p) 9.15
8.81 8.66 8.99
9.08 9.36 9.55
10.91 11.19 11.38
0.11 0.54 0.40
2.10 2.53 2.39
733.9 883.4 834.3
m[g]
EZPE[h]
s99.9%[j]
s[i] 213
2572.9 1286.5 1.3110 9.0610213 213 2572.9 1286.5 2.7010 1.8610212 2572.9 1286.5 2.1310213 1.4710212
[a] The relative Gibbs free energy of the C*T* base pair (DGCT 5 0; T 5 298.15 K), kcal mol21. [b] The relative electronic energy of the C*T* base pair (DECT 5 0), kcal mol21. [c] The Gibbs free energy of activation for the forward reaction of tautomerization (T 5 298.15 K), kcal mol21. [d] The activation electronic energy for the forward reaction of tautomerization, kcal mol21 [e] The Gibbs free energy of activation for the reverse reaction of tautomerization (T 5 298.15 K), kcal mol21. [f ] The activation electronic energy for the reverse reaction of tautomerization, kcal mol21. [g] The frequency of the vibrational mode in the C*T* base pair, which becomes imaginary in the TSCT $ C*T* of tautomerization obtained at the B3LYP/6-31111G(d,p) level of geometry optimization, cm21. [h] The zero-point vibrational energy associated with this normal mode, cm21. [i] The lifetime of the C*T* base pair, s. [j] The time necessary to reach 99.9% of the equilibrium concentration of the CT reactant and the C*T* product of the CT $ C*T* tautomerization reaction via the DPT, s.
Journal of Computational Chemistry 2013, DOI: 10.1002/jcc.23412
5
FULL PAPER
WWW.C-CHEM.ORG
Figure 1. Geometric structures of the nine key points describing the evolution of the CT $ C*T* tautomerization via the DPT along the IRC obtained at the B3LYP/6-31111G(d,p) level of theory in vacuo. Coordinates of each key point are presented below it. The dotted lines indicate AH…B H-bonds and O2…O2 vdW contact, whereas continuous lines show covalent bonds (their lengths are presented in angstroms). Carbon atoms are in light-blue, nitrogen in dark-blue, oxygen in red, and hydrogen in grey.
structure.[49,67,96] The time s99.9% necessary to reach 99.9% of the equilibrium concentration of the starting CT and the final C*T* base pair is equal to 1.47 3 10212 s (Table 3). The small value (
