1826
J. Phys. Chem. A 2010, 114, 1826–1834
Hydrolytic Deamination of 5,6-Dihydrocytosine in a Protic Medium: A Theoretical Study Vanessa Labet,*,† Christophe Morell,† Thierry Douki,† Jean Cadet,† Leif A. Eriksson,‡ and Andre´ Grand† Laboratoire “Le´sions des Acides Nucle´iques”, INAC/SCIBsUMR-E no. 3 CEA-UJF, CEA Grenoble, 17 rue ¨ rebro Life des Martyrs, F-38 054 Grenoble cedex 9, France, and School of Science and Technology and O ¨ ¨ Science Center, Orebro UniVersity, 701 82 Orebro, Sweden ReceiVed: May 26, 2009; ReVised Manuscript ReceiVed: NoVember 25, 2009
The mechanism for the deamination reaction of 5,6-dihydrocytosine with H2O in a protic medium was investigated by DFT calculations at the B3LYP/6-311G(d,p) level of theory as a model reaction for the deamination reaction of 5,6-saturated cytosine derivatives. Two pathways were found. Pathway dhA, which can explain the deamination in a protic medium at acidic pH, and pathway dhBt, more representative of the reaction in a protic medium at neutral pH. Pathway dhA is a two-step mechanism initiated by the nucleophilic addition of a water molecule to carbon C4 of N3-protonated 5,6-dihydrocytosine with the assistance of a second water molecule, followed by elimination of an ammonium cation to form 5,6-dihydrouracil. The nucleophilic addition is rate-determining with an activation free energy of 116.0 kJ/mol in aqueous solution. Pathway dhBt is a four-step mechanism which starts with the water-assisted tautomerization of 5,6dihydrocytosine to form the imino tautomer. This intermediate undergoes nucleophilic addition of water to carbon C4, which after protonation eliminates an ammonium cation, as in pathway dhA. The nucleophilic addition is again rate-determining, with an activation free energy of 113.3 kJ/mol in aqueous solution. The latter value is about 25 kJ/mol lower than its counterpart for cytosine, in agreement with the experimental observation that 5,6-saturated cytosine derivatives exhibit a much shorter lifetime in aqueous solution than their unsaturated counterparts. The evaluation of reactivity indices derived from conceptual DFT leads to the conclusion that this lower activation free energy can be attributed to a larger local electrophilic power of carbon C4 in 5,6-saturated derivatives. 1. Introduction Modification of the chemical structure of DNA is a major deleterious process in cells since it may lead to lethality and appearance of mutations. DNA damage is induced by a variety of physical (ionizing and UV radiation) and chemical (oxidative stress, chemicals) agents. In cells, spontaneous hydrolytic reactions are also at the origin of modifications in the chemical structure of DNA.1 For example, the N-glycosidic bond that links 2-deoxyribose units to the bases may be cleaved mostly at purine nucleotides, generating an abasic site. Another reaction involving hydrolytic processes is deamination of cytosine. The net result of this process is the substitution of the exocyclic amino group of cytosine by a hydroxyl function, leading to the enolic form of uracil. This spontaneous conversion of cytosine into uracil is quite frequent in cells as shown by the predominance of C to T transitions in the spontaneous mutation spectra. It is also worth noting that cells are equipped with uracil N-glycosylase, an enzyme that removes uracil from DNA. In spite of its undoubtful biological relevance, cytosine deamination is a very slow and inefficient process. In doublestranded DNA, the half-life of cytosine was estimated by a biochemical approach to be more than 30 centuries under physiological conditions.2 This value dropped by 2 orders of magnitude in single-stranded DNA. Data were also reported for monomers such as the free base and the ribo- and 2′* To whom correspondence should be addressed. Phone: +33 438 785 376. Fax: +33 438 785 090. E-mail:
[email protected]. † CEA Grenoble. ‡ ¨ Orebro University.
deoxyribonucleosides. These studies had to be conducted under conditions of extreme pH and temperature to accurately quantify the extent of deamination, albeit still remaining a slow process. For example, the rate constant for the deamination of cytidine is 3 × 10-6 s-1 at 37 °C in 0.01 M HCl3 and 10-6 s-1 at pH 6 and 90 °C.4 However, saturation of the 5,6-double bond of cytosine drastically increases the deamination rate. This was first observed for 5,6-dihydrocytosine, which exhibits a half-life of less than 3 h in aqueous solution.5 The enhanced deamination rate of 5,6-saturated cytosine derivatives was also observed in oxidation6-8 and hydration9,10 products. Deamination rate constants ranging between 5 × 10-5 and 10-7 s-1 under physiological conditions were similarly observed for saturated cytosine moieties in dimeric photoproducts produced upon UV irradiation, such as cyclobutane pyrimidine dimers and pyrimidine (6-4) pyrimidone photoproducts.11-15 Interestingly, deamination of UV-induced saturated dimeric cytosine photoproducts has been proposed16,17 to account for the typical C to T mutations observed in sunlight-induced skin tumors.18,19 A better understanding of the mechanism of deamination of 5,6-saturated cytosine derivatives is thus required to explain the bulk of these data. In the present work, emphasis was placed on the deamination of 5,6-dihydrocytosine that can be considered as a model for 5,6-saturated cytosine derivatives. We used quantum chemistry calculations to describe the reaction, from the initial attack of the water molecule to the final formation of uracil. On the basis of a previous work on unmodified cytosine20 where two pathways were proposed to explain the spontaneous deamination of cytosine in a protic medium (pathway A and
10.1021/jp9049044 2010 American Chemical Society Published on Web 01/07/2010
Hydrolytic Deamination of 5,6-Dihydrocytosine
J. Phys. Chem. A, Vol. 114, No. 4, 2010 1827
Figure 2. Atomic labeling used for the description of (a) pathway dhA and (b) pathways dhB and dhBt.
Figure 1. Reaction pathways involved in the spontaneous deamination of cytosine (red labels) and 5,6-dihydrocytosine (black labels) in a protic medium.
pathway B, Figure 1), three pathways (dhA, dhB, and dhBt) were studied in the case of 5,6-dihydrocytosine. They are displayed in Figure 1. Pathways dhA and dhB are the 5,6saturated counterparts of pathways A and B. In pathway dhA, the N3-protonated form of 5,6-dihydrocytosine undergoes a nucleophilic attack by a water molecule on the carbon holding the amino exocyclic group, with the assistance of a second water molecule. An ammoniac molecule is then released, forming a protonated form of 5,6-dihydrouracil (dhI2), in acidic equilibrium with an ammonium cation and 5,6-dihydrouracil (dhP1). In pathway dhB, 5,6-dihydrocytosine is directly attacked by the water dimer. A neutral tetrahedral intermediate is formed (dhI3) which is protonated to lead to intermediate dhI1, like in pathway dhA. The end of the pathway is the same as that previously proposed. Pathway dhBt is quite similar to pathway dhB, differing only by the formation of the tetrahedral intermediate dhI3. In pathway dhB, the latter intermediate is formed in a single step from dhR2, whereas in pathway dhBt there is a first tautomerization step whereafter the imine tautomer (dHR2t), an intermediate between 5,6-dihydrocytosine and its protonated form, undergoes the nucleophilic attack to form dhI3 via dHTS4t. The atomic labeling used in the following for the study of pathways dhA, dhB, and dhBt is indicated in Figure 2. 2. Computational Details All the calculations were performed using the Gaussian 03 package21 at the B3LYP/6-311G(d,p)22-26 level. 2.1. Computational Investigation of Reaction Pathways. As a first step, molecular systems of interest were geometry optimized in a vacuum. After full optimization, the different stationary points were characterized as minima or transition states by computing the vibrational frequencies within the
harmonic approximation at the same level of theory. Thermal data were extracted to obtain the different thermodynamic functions of reaction and the corresponding activation parameters at 298.15 K and 1 atm. The B3LYP functional is known to predict with moderate accuracy the barrier heights accompanying H-transfers. Thus, single-point MP2 calculations with the same basis set were also performed to check the barrier heights. It appeared that as expected B3LYP underestimates the barrier heights but that the general energy trend was not affected. Single-point energy calculations were performed on the optimized geometries including a polarized continuum model in the integral equation formalism (IEF-PCM, hereafter named PCM)27-29 with dielectric constant ε ) 78.39 to simulate the bulk effect of a polar environment. This value is clearly not appropriate to simulate the quite apolar environment inside duplex DNA, but it is adapted to reproduce that of nucleoside 5′-monophosphate in an aqueous solvent, for which was experimentally observed a difference of reactivity between cytosine and 1-methyl-5,6-dihydrocytosine.4,30 The cavity was built using the united atom topologic model applied to the atomic radii of the UFF force field, with an average area of 0.2 Å2 for the tesserae generated on each sphere. Moreover, the default cavity was modified by adding individual spheres to all hydrogen atoms linked to nitrogen and oxygen atoms, using the keyword SPHEREONH. Reactions paths in a vacuum were characterized at the same level of theory by performing intrinsic reaction coordinate (IRC) calculations from each gas-phase optimized transition structure to ensure that they connected the appropriate reactants and products. 2.2. Optimizations of the Reactive Complexes dhR1, dhR2, and dhR2t. The reactive forms of the three pathways, dhR1, dHR2, and dhR2t, are complexes involving a derivative of 5,6-dihydrocytosine in interaction with two water molecules. To get the correct relative orientation between the three molecules, dhR1 and dhR2t had to be optimized under constraint, as outlined in our recent paper31 dealing with spontaneous deamination of 5-methylcytosine in a protic medium. This is based on the observation of the reaction force32 profiles. For the details, the reader is referred to ref 30. The geometries of dhR1′ and dhR2t′ obtained with this method are represented in Figure 3, the prime indicating that the optimizations are performed under constraint, as well as that of dhR2. The three reactant complexes are stabilized by three H-bonds. One water molecule is in the plane defined by atoms N3, C4, and N4 and forms two H-bonds with 5,6dihydrocytosine. The second water molecule is above this plane and interacts with the first water molecule via one H-bond.
1828
J. Phys. Chem. A, Vol. 114, No. 4, 2010
Labet et al.
Figure 3. Optimized geometries of the dhR1′, dhR2, and dhR2t′ complexes in a vacuum. Distances are given in angstroms. dhR2 was optimized without imposing any constraint. For dhR1′, the distance between oxygen O8 and nitrogen N3 was fixed to 2.840 Å and the dihedral angle formed by oxygen O8, nitrogen N3, nitrogen N1, and carbon C6 to -82.3°. For dhR2t′, the distance between oxygen O8 and nitrogen N4 was fixed to 3.322 Å.
2.3. Evaluation of Reactivity Indices. Several reactivity indices derived from conceptual DFT33-38 and allowing a measure of electrophilicity of molecular systems39,38 were used in this work to rationalize the difference of reactivity between cytosine and 5,6-dihydrocytosine and their imine tautomers and N3-protonated forms toward hydrolytic deamination. They were calculated from geometries optimized on one hand in a vacuum and on the other hand in an aqueous solvent modelized as described in section 2.1. Some of the indices used in this work are global; others are local. Moreover, some are more accurate to describe electrostatic control (or charge control), whereas others are more adapted to describe an electron-transfer control (or frontier orbital control).40-43 2.3.1. Definition and Computation of the Global ReactiWity Indices Used in This Work. 2.3.1.1. Chemical Potential (µ). Within the conceptual DFT, the chemical potential µ has been defined as33,34,44
∂E µ) ∂N
( )
V(b) r
(1)
It measures the escape tendency of the electronic cloud from equilibrium.45 Consequently, a good electrophilic reactant must have a low chemical potential µ. 2.3.1.2. Absolute Hardness (η). The absolute hardness η, which can be viewed as the resistance of a molecular system to a charge transfer,46-48 is defined as33,34,45
η)
( ) 2
∂E ∂N2
(2) V(b) r
For a given µ, the smaller the absolute hardness, the bigger the electrophilicity of a molecular system.49
2.3.1.3. Global Electrophilicity Index (ω). The global electrophilicity index ω is defined as50,51
ω)
µ2 2η
(3)
By analogy with the equation of power in classical electricity (W ) V2/R, with W the power, V the potential difference, and R the resistance), ω can be considered as the measure of the electrophilic power of a molecular system. Assuming that the chemical potential is often negative, a large electrophilic power is consistent with a low (very negative) chemical potential and a small absolute hardness, as already mentioned. 2.3.1.4. Computational EValuation. µ, η, and ω were calculated using approximate expressions deriving from the use of a finite difference approximation and the Janak theorem:52,53
1 µ ≈ (εLUMO + εHOMO) 2
(4)
η ≈ εLUMO - εHOMO
(5)
ω≈
2 1 (εLUMO + εHOMO) 8 εLUMO - εHOMO
(6)
where εHOMO and εLUMO are the energies of the highest occupied (HOMO) and the lowest unoccupied (LUMO) molecular orbitals, respectively. 2.3.2. Definition and Computation of the Local ReactiWity Indices Used in This Work. All the previous indices provide information on the reactivity of a molecular system as a whole. However, they are not relevant for assessing differences in
Hydrolytic Deamination of 5,6-Dihydrocytosine
J. Phys. Chem. A, Vol. 114, No. 4, 2010 1829
reactivity within the system. For this, local indices were used: the net atomic charges to describe reactivity under electrostatic control and three other indices more accurate to describe reactivity under electron-transfer control. 2.3.2.1. Net Atomic Charges. The net atomic charges obtained from the electrostatic potential provide a good description of the way in which the electron density distribution of a molecular system interacts with other molecules. b)). The local elec2.3.2.2. Local Electrophilicity Index (ω+(r b), which is a local version of the global trophilicity index ω+(r electrophilicity index ω, is defined as54
ω+(b) r ) ωf+(b) r
(7)
where f+(r b) is the Fukui function used when the system undergoes a nucleophilic attack:55,56
f+(b) r )
[ ∂F(∂Nb)r ]
+
(8)
V
with F the electronic density and N the number of electrons in b) measures the intramolecular reactivthe molecular system. f+(r ity at site b r toward a nucleophilic reagent. The bigger the f(r b), the more reactive the site with respect to other sites of the same molecule toward a nucleophilic attack. Since ω measures the b) takes into account both global electrophilicity power, ω+(r global electrophilicity and local intramolecular reactivity. 2.3.2.3. Excess Electrophilicity (∆ω(r b)). The excess electrophilicity ∆ω(r b) is defined as57
∆ω(b) r ) ω∆f(b) r
(9)
where ∆f+(r b) is the dual descriptor,58,59 a third-order response 60 function, defined as
∆f(b) r )
( ) ∂2F(b) r ∂N2
(10)
V
If ∆f(r b) > 0, the site is more electrophilic than nucleophilic and the reverse if ∆f(r b) < 0. ∆ω(r b), like ω+(r b), takes into account both global electrophilicity and local intramolecular reactivity and is consequently a good intermolecular descriptor. Since ω is always positive, the more positive the ∆ω(r b), the more favored the site to undergo nucleophilic attack, whereas the more negative the ∆ω(r b), the more favored the site to undergo electrophilic attack. 2.3.2.4. Another Local ReactiVity Index DeriVed from the EquiValent of the Dual Descriptor in the Grand Canonical Ensemble. Another reactivity indicator based on the dual descriptor ∆f(r b) was introduced in a recent paper61 to allow information to be gained on the relative reactivity of molecules with different sizes. This indicator is the equivalent of the ∆f(r b) function in the grand canonical ensemble and can be expressed as62
s(2)(b) r )
{( ) }
∆f(b) r f(b) r ∂η - 3 2 ∂N η η
V(b) r
(11)
This descriptor, contrary to ∆ω(r b) and ω+(r b) is size extensive. Consequently, it scales correctly the philicity with respect to system size and can be used as an appropriate tool to compare
the relative reactivity of cytosine and 5,6-dihydrocytosine.31 The b), the more favored a nucleophilic attack more positive the s(2)(r on that site. Assuming that the second term of the right-hand side of eq 11 is negligible with respect to the first one, the values of ∆f(r b)/η2 were used in our study to compare the electrophilicity of cytosine and 5,6-dihydrocytosine and that of their N3protonated forms. 2.3.2.5. Computational EValuation. The partial charges on the nuclei have been computed with the ChelpG63 method. b), ∆ω(r b), and ∆f(r b)/η2 were evaluated through their ω+(r condensed-to-atom version where the atom of interest is carbon C4. To calculate these, the following formulas were used: + + ωC4 ) ωfC4
(12)
∆ωC4 ) ω∆fC4
(13)
+ fC4 ) pC4(N + 1) - pC4(N)
(14)
∆fC4 ) pC4(N + 1) + pC4(N - 1) - 2pC4(N)
(15)
where pC4(N + 1), pC4(N - 1), and pC4(N) are, respectively, the electronic populations of carbon C4 when an electron is added to the molecular system, when an electron is removed, and in the normal system. 3. Results and Discussion 3.1. Energetics of Pathways dhA, dhB, and dhBt. Pathways dhA and dhB (Figure 1) have been investigated to assess how the loss of aromaticity due to the C5-C6 saturation may affect the way 5,6-dihydrocytosine undergoes deamination with respect to cytosine. All stationary points involved in pathway dhA were optimized, except dhI2, which could not be located. If present, dhI2 would be a complex among a water molecule, an ammoniac molecule, and a protonated form of 5,6-dihydrouracil. The instability of such systems has been described previously in the literature.31,64 It can be assumed that the level of calculation is not good enough to characterize this critical point. Since it appears quite clearly that dhI2 does not play a major role either in the reaction mechanism or in the difference of reactivity between cytosine and 5,6-dihydrocytosine, this point was not studied further. Consequently, pathway dhA is a two-step mechanism which involves initially a nucleophilic attack of a water molecule on carbon C4 of N3-protonated 5,6-dihydrocytosine, with the assistance of a second water molecule. This gives dhI1, a positively charged complex between a water molecule and a protonated 5,6-dihydrocytosine derivative in which carbon C4 is tetrahedral. An ammonium cation is then released from the 5,6-dihydrocytosine compound to immediately give dhP1, an H-bonded complex among this ion, a water molecule, and 5,6-dihydrouracil. The study of pathway dhB leads to the conclusion that it is not a valid pathway for this system. All attempts aimed at finding a transition-state structure linking dhR2 and dhI3 failed. However, it was possible to find a transition state linking dhI3 and a tautomeric form of dhR2. This suggests that, in the case of 5,6-dihydrocytosine, pathway B, present in cytosine and other 5,6-unsaturated cytosine derivatives such as 5-methylcytosine,21 is replaced by pathway dhBt, a four-step mechanism with an initial tautomerization step leading to the imine tautomer of 5,6dihydrocytosine in interaction with a water dimer (dhR2t′). A
1830
J. Phys. Chem. A, Vol. 114, No. 4, 2010
Labet et al.
TABLE 1: Relative Energies (kJ/mol) in a Vacuum, Free Energies (kJ/mol) in a Vacuum, and Relative Free Energies (kJ/mol) in an Aqueous Solvent for the Stationary Points of Pathways dhA and dhBt of 5,6-Dihydrocytosine and Pathways A and B of Cytosine31 pathway dhA
a
dhBta
a
system dhR1′ dhTS1 dhI1 dhTS2 dhI2 dhP1 dhR2 + H+ dhTSt + H+ dhR2t′ + H+ dhTS4 + H+ dhI3 + H+ dhI1 dhTS2 dhI2 dhP1
∆Evac
∆Gvac
∆Gaq
0.0 138.1 76.8 84.1
0.0 154.8 89.4 95.3
0.0 116.0 70.7 93.9
-55.6 0.0 42.6 -11.7 100.4 14.6 -926.6 -919.3
-48.6 0.0 51.7 -1.4 118.0 27.1 -890.1 -884.1
-62.5 0.0 44.0 3.4 113.3 28.1 22.7 45.9
-1059.0
-1028.1
-110.5
pathway b
A
Bb
system
∆Evac
∆Gvac
∆Gaq
R1′ TS1 I1 TS2 I2 P1 R2′ + H+
0.0 145.1 77.9 84.9 36.8 -43.2 0.0
0.0 163.8 87.9 93.9 20.0 -47.8 0.0
0.0 136.0 64.2 82.3 35.0 -65.2 0.0
TS4 + H+ I3 + H+ I1 TS2 I2 P1
121.4 62.1 -908.7 -901.6 -949.7 -1029.7
142.9 74.9 -865.6 -859.6 -933.5 -1001.3
138.5 79.3 44.0 62.1 14.8 -85.4
b
5,6-Dihydrocytosine, this work. Cytosine, ref 31.
water molecule attacks carbon C4 of the 5,6-dihydrocytosine tautomer with the assistance of the second water molecule. This results in the formation of a neutral 5,6-dihydrocytosine derivative in which C4 is tetrahedral, in interaction with a water molecule (dhI3). The cytosine derivative is then protonated on nitrogen N4 to allow the release of an ammonium cation in the final step. This gives dhP1, the same product as in pathway dhA. The relative energies in a vacuum and the relative Gibbs free energies in a vacuum and in an aqueous solvent for the stationary points of pathways dhA and dhBt are reported in Table 1. In addition, the values obtained in the case of unmodified cytosine, i.e., pathways A and B, are provided. At this point, it is worth noting that the authors are aware of the fact that the B3LYP functional with which were calculated those energies is known to underestimate by 20 kJ/mol the barrier heights of reactions involving H-transfers, as some of those involved in the studied pathways. This was checked by computing single-point MP2 calculations on B3LYP-optimized geometries. Since we are more interested in understanding reactivity differences than in computing reliable energy barriers, the rest of the discussion will be led with respect to the B3LYP energies. Nonetheless, the reader is invited to consider the absolute energy barriers with caution. In pathways dhA and dhBt, the nucleophilic addition of water is the rate-determining step, as in pathways A and B. It appears that the C5-C6 saturation facilitates these steps, a fact which is in agreement with the experimental observations of faster deamination rates mentioned in the Introduction. Indeed, the activation free energy for the dhR1′f dhI1 step is 9 kJ/mol lower than that for R1′f I1 in a vacuum and 20.0 kJ/mol in aqueous solution. As for the rate-determining step of pathway dhBt, this is associated with a 24.9 kJ/mol lower free energy barrier in a vacuum and 25.2 kJ/mol in an aqueous solvent. Such differences in activation free energies correspond to a ratio on the order of 3 × 105 between the rate constants, according to the Eyring equation.65 Although this factor is about 1 order of magnitude lower than that suggested from experimental studies,4,30 pathways dhA and dhBt do indeed explain why 5,6saturated cytosine derivatives have half-lives much shorter than that of unmodified cytosine. As noted for unmodified cytosine, pathway dhA is much easier in an aqueous solvent than in a vacuum. This can be explained by the fact that dhTS1 is more polar than dhR1′ (see
TABLE 2: Dipole Moments (µD, D) of the Stationary Points Involved in Pathways dhA, dhBt, A, and B in a Vacuum
a
systema
µD
systemb
µD
dhR1′ dhTS1 dhI1 dhTS2
6.4 8.9 4.6 3.7
dhP1 dhR2 dhTSt dhR2t′ dhTS4 dhI3
6.0 4.0 6.9 4.9 5.3 4.3
R1′ TS1 I1 TS2 I2 P1 R2′
5.2 7.4 6.4 6.0 3.6 11.3 4.8
TS4 I3
4.9 4.0
5,6-Dihydrocytosine, this work. b Cytosine, ref 31.
Figure 4. Relative free energies (kJ/mol) along pathways dhA and dhBt in aqueous solution.
Table 2). On the contrary, pathway dhBt is slightly more difficult in an aqueous solution than in a vacuum. The consequence is that in an aqueous solvent, pathway dhBt appears somewhat favored over pathway dhA, as seen in Table 1 and Figure 4. This result is in contradiction with experimental results obtained by Slae and Shapiro 30 years ago,26 who observed that, at 20 °C, 1-methyl-5,6-dihydrocytosine deaminates about 1.6 ( 0.1-fold faster at pH 0.4 than at pH 8.55-8.95. It should be noted, however, that according to the values reported in Table 1, the free energy barrier is almost the same in aqueous solution for pathways dhA and dhBt (difference of only 2.7 kJ/mol).
Hydrolytic Deamination of 5,6-Dihydrocytosine
J. Phys. Chem. A, Vol. 114, No. 4, 2010 1831
Figure 6. Optimized geometry of dhTSt. Distances are expressed in angstroms. The blue arrows are indicative of the normal mode of the imaginary frequency of the transition state.
Figure 5. Optimized geometries of dhTS1, TS1, dhI1, and I1 in a vacuum. Distances are in angstroms. The blue arrows on dhTS1 are indicative of the normal mode of the imaginary frequency of the transition state.
Also worth noting is that, in aqueous solution, pathway dhBt is 25.1 kJ/mol more exergonic than pathway B. Hence, not only is the rate-determining reaction easier for 5,6-dihydrocytosine with respect to unmodified cytosine, but the driving force of the whole reaction is also larger. 3.2. Geometries of the Stationary Points. After the previous subsection dealing with the comparison of the pathway energetics for 5,6-dihydrocytosine and unmodified cytosine, this subsection is devoted to the comparison of the geometries of the different stationary points involved. 3.2.1. Comparison of Transition States and Products of Pathway A and Pathway dhA Rate-Determining Steps. The first step of pathway dhA implies a concerted mechanism with the creation of a C-O single bond between carbon C4 of the N3-protonated form of 5,6-dihydrocytosine and O8 of one water molecule, a proton transfer from O8 to O7, and another proton transfer from O7 to N4. This results in the formation of a positively charged derivative of 5,6-dihydrocytosine in which C4 is tetrahedral (dhI1, Figure 5). The transition state dhTS1, whose geometry is displayed in Figure 5, is a six-membered transition state (C4-O8-H8a-O7-H7a-N4). Its structure is very similar to that of its counterpart in the case of pathway A (TS1, Figure 5), with differences in the distances being less than 0.02 Å. At the transition state, in the case of N3-protonated 5,6-dihydrocytosine as well as for N3-protonated cytosine, the proton transfer of H8a is much more advanced than that of H7a, indicating that the concerted mechanism is quite asynchronous. The intermediate dhI1 is, like I1, a hydrogen-bonded complex between the positively charged tetrahedral derivative of (5,6dihydro)cytosine and a water molecule. Because of the higher flexibility of the C5-C6 saturated derivative of cytosine, the water molecule can involve all three atoms in hydrogen bonds in dhI1, whereas I1 is stabilized by only two H-bonds. 3.2.2. Transition State and Product of the Tautomerization Step of Pathway dhBt. The transition state of the first step of pathway dhB, a tautomerization step which has no equivalence in pathway B, is displayed in Figure 6. Albeit in dhR2 5,6-
Figure 7. Optimized geometries of dhTS4, TS4, dhI3, and I3. Distances are given in angstroms. The blue arrows on dhTS4 are indicative of the normal mode of the imaginary frequency of the transition state.
dihydrocytosine interacts with two water molecules, only one participates actively in the tautomerization. Consequently, this step corresponds to a double proton transfer: H7a is transferred from O7 to N3 and H4a from N4 to O7. dhTSt is a sixmembered transition state. The distances involved indicate that the two proton transfers are quite synchronous. It can be noted that although the second water molecule is not active in this step, the H-bond between the two water molecules at the transition state is stronger than in dhR2 and dhR2t′ since the distance between O7 and H8a is shorter in dhTSt (1.788 Å) than in dhR2 (1.931 Å) and dhR2t′ (1.973 Å). This can be explained by the fact that during tautomerization the electronic density on O7 is increased, which hence attracts the weakly positive hydrogen. 3.2.3. Comparison of Transition States and Products of Pathway B and Pathway dhBt Rate-Determining Steps. The second step of pathway dhB is similar to the first step of pathway B. It involves a first proton transfer from a water molecule to a nitrogen atom of the cytosine derivative (H4a transferred from O7 to N4 in the case of 5,6-dihydrocytosine, H7a transferred from O7 to N3 in the case of cytosine), a second proton transfer (H8a) between the two water molecules, and the creation of a covalent bond between carbon C4 and oxygen O8. This leads
1832
J. Phys. Chem. A, Vol. 114, No. 4, 2010
Labet et al. TABLE 3: Kohn-Sham Frontier Orbital Energies (εLUMO and εHOMO, eV), Chemical Potentials (µ, eV), Absolute Hardnesses (η, eV), and Global Electrophilicity Indices (ω, eV) in a Vacuum (a) and in an Aqueous Solvent (b) molecular system
εLUMO
εHOMO
µ
η
ω
(a) Vacuum cytosine -1.03 5,6-dihydrocytosine -0.46 cytosine imine tautomer -1.09 5,6-dihydrocytosine imine tautomer -0.16 N3-protonated cytosine -6.64 N3-protonated 5,6-dihydrocytosine -6.20
-6.39 -6.50 -6.50 -7.02 -11.86 -12.05
-3.71 -3.48 -3.80 -3.59 -9.25 -9.13
5.36 6.04 5.41 6.86 5.22 5.85
1.29 1.00 1.33 0.94 8.19 7.12
-3.74 -3.70 -3.62 -3.66 -4.86 -4.84
5.58 6.26 5.50 6.88 5.36 6.04
1.25 1.09 1.19 0.97 2.20 1.94
(b) Aqueous Solvent cytosine -0.95 -6.52 5,6-dihydrocytosine -0.57 -6.83 cytosine imine tautomer -0.87 -6.37 5,6-dihydrocytosine imine tautomer -0.22 -7.10 N3-protonated cytosine -2.18 -7.53 N3-protonated 5,6-dihydrocytosine -1.82 -7.86
Figure 8. Optimized geometries of dhTS2, TS2, dhP1, and P1. Distances are given in angstroms. The blue arrows on dhTS2 are indicative of the normal mode of the imaginary frequency of the transition state.
to the formation of a neutral molecular system, dhI3, whose geometry is displayed in Figure 7, in which C4 is tetrahedral. The transition state associated with this elementary step, dhTS4 (Figure 7), has a structure close to that of TS4, in which all the atoms involved in the step are organized in three six-membered rings: C4-N3-H7a-O7-H4a-N4, C4-N3-H7a-O7-H8aO8, and C4-N4-H4a-O7-H8a-O8. Due to the fact that, in pathway dhBt, H4a is transferred from O7 to N4 whereas, in pathway B, H7a is transferred from O7 to N3, the distances between O7 and H7a on one hand and O7 and H4a on the other hand are somewhat different in dhTS4 and TS4 (in dhTS4, dO7-H4a ) 1.795 Å, dO7-H7a ) 1.857 Å; in TS4, dO7-H4a ) 1.881 Å, dO7-H7a ) 1.783 Å). The remaining distances are more similar and differ by less than 0.02 Å. Like the first step of pathway dhA, the second step of pathway dhBt appears quite asynchronous since the distance between H4a and O7 in dhTS4 is much longer than that between H8a and O8 (in dhTS4, dO7-H4a ) 1.795 Å, dO8H8a ) 1.261 Å). The H8a proton transfer and the C-O bond creation occur later in dhTS4 with respect to TS4 since H8a is closer to O8 in dhTS4 than in TS4. This is also supported by the fact that O8 is further from C4 in dhTS4 than in TS4 (in dhTS4, dO8-H8a ) 1.261 Å, dO8-C4 ) 2.218 Å; in TS4, dO8-H8a ) 1.283 Å, dO8-C4 ) 2.204 Å). It may be noted that dhI3 and I3 have very similar structures. The main difference concerns the distance between the oxygen of the water molecule and the two hydrogens (H7a and H4a) of the cytosine derivative involved in H-bonds with it (in dhI3, dO7-H7a ) 2.215 Å, dO7-H4a ) 2.290 Å; in I3, dO7-H7a ) 2.159 Å, dO7-H4a ) 2.268 Å). 3.2.4. Transition States and Products InWolWed in NH3 Release. In dhTS2, whose geometry is displayed in Figure 8, as well as that of TS2, dhP1, and P1, the distance between C4 and N4 is 1.990 Å. This is about 0.13 Å longer than in TS2, whereas in I1 and dhI1 they are almost the same: 1.574 Å in dhI1 and 1.584 Å in I1. It may be added that dhTS2 is a later transition state than TS2. According to the Hammond postulate,
this is concordant with the fact that the C4-N4 bond breakage is a little more difficult for 5,6-dihydrocytosine than for unmodified cytosine (cf. Table 1). It can be noted that, in dhTS2, the water molecule is involved in three H-bonds, whereas only two are implicated in the stabilization of TS2. This is due to a slightly different orientation of the N4 group resulting from the C5-C6 saturation. In the same way, dhP1 shows one additional H-bond with respect to P1, between the water molecule and the ammonium cation. 3.3. Comparison of the Electronic Properties of Cytosine and 5,6-Dihydrocytosine, Including Those of Their Imine Tautomers and N3-Protonated Forms. To understand why pathways dhA and dhBt are easier than pathways A and B, respectively, and why pathway dhB does not exist, this section is devoted to the comparison of the electronic properties of cytosine and 5,6-dihydrocytosine together with those of their imine tautomers and N3-protonated forms. Since the deamination reaction involves electrophilic properties of the cytosine derivative, it was decided to focus on global and local reactivity indices from conceptual DFT, which measure these properties. 3.3.1. Global ReactiWity Indices. In Table 3 we report the energies of the frontier orbitals of cytosine and 5,6-dihydrocytosine and their imine tautomers and N3-protonated forms as well as their chemical potentials, absolute hardnesses, and global electrophilicity indices. The C5-C6 bond saturation stabilizes the HOMO and destabilizes the LUMO. This holds true both in a vacuum and in an aqueous solution for the canonical cytosine and its imine tautomer and N3-protonated form. As a consequence, 5,6dihydrocytosine and its imine tautomer and N3-protonated 5,6dihydrocytosine are harder than their cytosine counterparts by 0.6-0.7, 1.4-1.5, and 0.6-0.7 eV, respectively. They are more resistant to electronic transfer, which is in agreement with the fact that in general saturated molecules are less reactive than unsaturated ones. In a vacuum, the chemical potential is less negative by 0.1-0.2 eV when the C5-C6 bond is saturated, whereas the effect is essentially negligible in an aqueous solution. This means that pyrimidine bases with a C5-C6 saturated bond are less electronegative than their unsaturated counterparts. Indeed, it is well-known that the more unsaturated a C-C bond, the more electronegative it is. The influence of the C5-C6 bond saturation on both chemical potential and absolute hardness implies that saturated derivatives have a smaller global electrophilicity index than the unsaturated ones. This can be easily understood since ω represents the energetic stabilization of a system that acquires enough electrons to be saturated and 5,6-dihydrocytosine is more saturated in electrons
Hydrolytic Deamination of 5,6-Dihydrocytosine
J. Phys. Chem. A, Vol. 114, No. 4, 2010 1833
TABLE 4: Electrostatic Charges on Carbon C4 (qC4), Local + Electrophilicity (ωC4 , eV), Excess Electrophilicity (∆ωC4, eV), and Grand Canonical Dual Descriptor (∆fC4/η2, 10-2 eV-2) Condensed to Carbon C4 in a Vacuum (a) and in an Aqueous Solvent (b) + ωC4
∆ωC4
∆fC4/η2
(a) Vacuum cytosine 0.91 5,6-dihydrocytosine 0.86 cytosine imine tautomer 0.73 5,6-dihydrocytosine imine tautomer 0.65 N3-protonated cytosine 0.85 N3-protonated 5,6-dihydrocytosine 0.75
0.27 0.40 0.11 0.26 2.57 3.99
0.27 0.34 0.13 0.30 2.64 3.69
0.71 0.92 0.33 0.67 1.18 1.51
(b) Aqueous cytosine 5,6-dihydrocytosine cytosine imine tautomer 5,6-dihydrocytosine imine tautomer N3-protonated cytosine N3-protonated 5,6-dihydrocytosine
0.34 0.59 0.21 0.49 0.72 1.22
0.32 0.58 0.18 0.49 0.67 1.13
0.81 1.34 0.49 1.05 1.06 1.59
molecular system
qC4
Solvent 1.07 0.93 0.81 0.73 0.90 0.77
than cytosine. The effect is more pronounced in a vacuum than in an aqueous solvent. As a general trend it may be concluded that, in their wholeness, 5,6-saturated derivatives of cytosine are less electrophilic than their unsaturated counterparts because of a larger absolute hardness and a higher (less negative) chemical potential. Since pathways dhA and dhBt appear easier than pathways A and B, respectively, we furthermore investigate the influence of the C5-C6 saturation not only on the global electrophilicity but also on the electrophilicity of carbon C4. 3.3.2. Local ReactiWity Indices. In Table 4 we report the atomic charge held by C4 in a vacuum and in an aqueous solvent for cytosine and 5,6-dihydrocytosine and their imine tautomers and N3-protonated forms as an indicator to describe charge control. The values of other reactivity descriptors derived from conceptual DFT adapted to describe frontier orbital control, such + , ∆ωC4, and ∆fC4/η2, are also reported in Table 4. It can as ωC4 be observed that C5-C6 saturation induces a slightly reduced positive charge held by carbon C4, by about 0.05-0.14 e, both in a vacuum and in an aqueous solvent. On the other hand, local + , ∆ωC4, and ∆fC4/η2 indicate that C4 is reactivity indices ωC4 more electrophilic in the 5,6-dihydrocytosine derivatives than in their counterparts derived from unmodified cytosine, even though these molecular systems are less electrophilic as a whole. Since pathways dhA and dhBt are easier than pathways A and B, respectively, it can be suspected that the deamination reaction is under frontier orbital control. Interestingly, the C4 carbons of the cytosine and 5,6dihydrocytosine imine tautomers appear much less electrophilic than those of cytosine and 5,6-dihydrocytosine, respectively. Consequently, the reason why pathway dhB is replaced by pathway dhBt is not related to the electrophilicity of C4. A hypothesis is that it is in relation with the proton transfers which take place during the rate-determining step. According to pathway dhB, N3 has to be protonated, whereas this holds for N4 in pathway dhBt. In Figure 9 we display the dual descriptor ∆f(r b) for 5,6dihydrocytosine (Figure 9a) and its imine tautomer (Figure 9b). Areas in red indicate that ∆f(r b) is positive and that the site hence is more electrophilic than nucleophilic; the reverse holds for green. It can be assumed that atoms able to be protonated correspond to nucleophilic sites and thus negative values of ∆f(r b). It appears from Figure 9a,b that nitrogen N3 of 5,6dihydrocytosine, contrary to nitrogen N4 of its imine tautomer,
Figure 9. ∆f(r b) calculated at the B3LYP/6-311G(d,p) level of theory for (a) 5,6-dihydrocytosine, (b) 5,6-dihydrocytosine imine tautomer, (c) N3-protonated 5,6-dihydrocytosine, (d) cytosine, (e) cytosine imine tautomer, and (f) N3-protonated cytosine. The red areas correspond to positive values of ∆f(r b) and the green areas to negative values.
is not favorable toward protonation. It is thus the imine tautomer which undergoes nucleophilic addition, even if the C4 carbon is less electrophilic than that of 5,6-dihydrocytosine. It is interesting to note that the situation is not the same for unmodified cytosine (cf. parts d and e of Figure 9); instead both N3 of cytosine and N4 of its imine tautomer can be protonated. Since the C4 atom of cytosine has a larger electrophilic power than its imine tautomer, it is the amino tautomer which is attacked by the water dimer. In parts c and f of Figure 9, which represent the dual descriptor for N3-protonated forms of cytosine and 5,6-dihydrocytosine, it appears that N4 is not favorable to protonation. Since the system is already positively charged, this is quite understandable and can be a reason why in a vacuum pathways dhA and A are more difficult compared to pathways dhBt and B. In our recent study,31 it had been assumed that this could arise from the fact that in a vacuum the chemical potential of N3-protonated cytosine is very different from that of a water molecule. Quite possibly, both these aspects have to be taken into account for a complete description. 4. Conclusions In this study, three pathways have been investigated to explain the spontaneous deamination of 5,6-dihydrocytosine in a protic medium. They are based upon pathways proposed in the literature for the spontaneous deamination of unmodified cytosine in the same type of environment. Two of the pathways, dhA and dhBt, display energy barriers able to explain why 5,6-saturated cytosine derivatives deaminate faster than unmodified cytosine. In both pathways dhA and dhBt the rate-determining step is the nucleophilic addition of a water molecule to carbon C4 with the assistance of a second water molecule. In pathway dhA, the electrophilic agent is N3protonated 5,6-dihydrocytosine, whereas in pathway dhBt, it is the 5,6-dihydrocytosine imine tautomer. The evaluation of global and local electrophilicity indices confirms that although 5,6-dihydrocytosine derivatives are less electrophilic than cytosine derivatives as a whole, their C4 carbon is more electrophilic. The third pathway, dhB, appears to be unlikely. It seems that 5,6-dihydrocytosine cannot undergo the addition of a water molecule because its nitrogen N3 is not sufficiently nucleophilic to accommodate the charge of a proton. This is the main difference between the way 5,6-saturated and 5,6-unsaturated neutral cytosine derivatives deaminate in a protic medium. According to our calculations, in a vacuum, pathway dhBt is easier than pathway dhA. This may be explained by the fact that N3-protonated 5,6-dihydrocytosine is
1834
J. Phys. Chem. A, Vol. 114, No. 4, 2010
less likely than neutral 5,6-dihydrocytosine to include the charge of one of the hydrogen atoms of the added water molecule. In an aqueous solvent, both pathways are associated with comparable energetic barriers (∆∆Gq ) 2.7 kJ/mol), with pathway dhBt being a little more favorable, which is in contradiction with experimental results. One explanation could be that the aqueous environment is not described well enough. Acknowledgment. L.A.E. thanks the Swedish Research Council (VR) and the Faculty of Science and Technology at ¨ rebro University for financial support. V.L., C.M., and A.G. O thank the Centre de Calcul Recherche et Technologie (CCRT) for allocation of computer time. Supporting Information Available: Molecular modeling coordinates of all stationary points involved in pathways dhA, dhBt, A, and B. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Lindahl, T. Nature 1993, 362, 709. (2) Frederico, L. A.; Kunkel, T. A.; Ramsay Shaw, B. Biochemistry 1990, 29, 2532. (3) Kusmierek, J.; Ka¨ppi, R.; Neuvonen, K.; Shugar, D.; Lo¨nnberg, H. Acta Chem. Scand. 1989, 43, 196. (4) Shapiro, R.; Klein, R. S. Biochemistry 1966, 5, 2358. (5) Green, M.; Cohen, S. S. J. Biol. Chem. 1957, 228, 601. (6) Decarroz, C.; Wagner, J. R.; van Lier, J. E.; Krishna, M. C.; Riesz, P.; Cadet, J. Int. J. Radiat. Biol. 1986, 50, 491. (7) Tremblay, S.; Douki, T.; Cadet, J.; Wagner, R. J. Biol. Chem. 1999, 274, 20833. (8) Wagner, J. R.; Hu, C. C.; Ames, B. N. Proc. Natl. Acad. Sci. U.S.A. 1992, 89, 3380. (9) Boorstein, R. J.; Hilbert, T. P.; Cunningham, R. P.; Teebor, G. W. Biochemistry 1990, 29, 10455. (10) Douki, T.; Vadesne-Bauer, G.; Cadet, J. Photochem. Photobiol. Sci. 2002, 1, 565. (11) Freeman, K. B.; Hariharan, P. V.; Johns, H. E. J. Mol. Biol. 1965, 13, 833. (12) Lemaire, D. G. E.; Ruzsicska, B. P. Photochem. Photobiol. 1993, 57, 755. (13) Lemaire, D. G. E.; Ruzsicska, B. P. Biochemistry 1993, 32, 2525. (14) Douki, T.; Cadet, J. J. Photochem. Photobiol., B 1992, 15, 199. (15) Liu, F.-T.; Yang, N. C. Biochemistry 1978, 17, 4865. (16) Peng, W.; Shaw, B. R. Biochemistry 1996, 35, 10172. (17) Jiang, N.; Taylor, J.-S. Biochemistry 1993, 32, 472. (18) Ziegler, A.; Leffel, D. J.; Kunala, S.; Sharma, H. W.; Shapiro, P. E.; Bale, A. E.; Brash, D. E. Proc. Natl. Acad. Sci. U.S.A. 1993, 90, 4216. (19) Brash, D. E.; Rudolph, J. A.; Simon, J. A.; Lin, A.; McKenna, G. J.; Baden, H. P.; Halperin, A. J.; Ponten, J Proc. Natl. Acad. Sci. U.S.A. 1991, 88, 10124. (20) Labet, V.; Grand, A.; Morell, C.; Cadet, J.; Eriksson, L. A. Theor. Chem. Acc. 2008, 120, 429. (21) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskortz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.;
Labet et al. Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chem, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision C.02; Gaussian, Inc.: Wallingford, CT, 2004. (22) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (23) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (24) McLean, A. D.; Chandler, G. S. J. Chem. Phys. 1980, 72, 5639. (25) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. J. Chem. Phys. 1980, 72, 650. (26) Frisch, M. J.; Pople, J. A.; Binkley, J. S. J. Chem. Phys. 1984, 80, 3265. (27) Cance`s, M. T.; Mennucci, B.; Tomasi, J. J. Chem. Phys. 1997, 107, 3032. (28) Mennucci, B.; Tomasi, J. J. Chem. Phys. 1997, 106, 5151. (29) Cossi, M.; Scalmani, G.; Rega, N.; Barone, V. J. Chem. Phys. 2002, 117, 43. (30) Slae, S.; Shapiro, R. J. Org. Chem. 1978, 43, 1721. (31) Labet, V.; Morell, C.; Cadet, J.; Eriksson, L. A.; Grand, A. J. Phys. Chem. A 2009, 113, 2524. (32) Toro-Labbe´, A. J. Phys. Chem. A 1999, 103, 4398. (33) Parr, R. G.; Yang, W. Density Functional Theory of Atoms and Molecules; Oxford University Press: New York, 1989. (34) Geerlings, P.; De Proft, F.; Langenaeker, W. Chem. ReV. 2003, 103, 1793. (35) Chermette, H. J. Comput. Chem. 1999, 20, 129. (36) Liu, S. B. Acta Phys.-Chim. Sin. 2009, 25, 590. (37) Gazquez, J. L. J. Mex. Chem. Soc. 2008, 52, 3. (38) Ayers, P. W.; Anderson, J. S. M.; Bartolotti, L. J. Int. J. Quantum Chem. 2005, 101, 520. (39) Chattaraj, P. K.; Sarkar, U.; Roy, D. R. Chem. ReV. 2006, 106, 2065. (40) Klopman, G. J. Am. Chem. Soc. 1968, 90, 223. (41) Melin, J.; Aparicio, F.; Subramanian, V.; Galvan, M.; Chattaraj, P. K. J. Phys. Chem. 2004, 108, 2487. (42) Ayers, P. W.; Parr, R. G.; Pearson, R. G. J. Chem. Phys. 2006, 124, 194107. (43) Anderson, J. S. M.; Melin, J.; Ayers, P. W. J. Chem. Theory Comput. 2007, 3, 358. (44) Parr, R. G.; Yang, W. Annu. ReV. Phys. Chem. 1995, 46, 701. (45) Parr, R. G.; Donnelly, R. A.; Levy, M.; Palke, W. E. J. Chem. Phys. 1978, 68, 3801. (46) Parr, R. G.; Pearson, R. G. J. Am. Chem. Soc. 1983, 105, 7512. (47) Pearson, R. G. J. Am. Chem. Soc. 1985, 107, 6801. (48) Pearson, R. G. J. Chem. Educ. 1987, 64, 561. (49) Ayers, P. W. Faraday Discuss. 2007, 135, 161. (50) Maynard, A. T.; Huang, M.; Rice, W. G.; Covel, D. G. Proc. Natl. Acad. Sci. U.S.A. 1998, 95, 11578. (51) Parr, R. G.; Szentpaˇly, L. V.; Liu, S. J. Am. Chem. Soc. 1999, 121, 1922. (52) Janak, J. F. Phys. ReV. B 1978, 18, 7165. (53) Cohen, A. J.; Mori-Sanchez, P.; Yang, W. T. Phys. ReV. B 2008, 77, 115123. (54) Chattaraj, P. K.; Maiti, B.; Sarkar, U. J. Phys. Chem. A 2003, 107, 4973. (55) Parr, R. G.; Yang, W. J. Am. Chem. Soc. 1984, 106, 4049. (56) Ayers, P. W.; Levy, M. Theor. Chem. Acc. 2000, 103, 353. (57) Padmanabhan, J.; Parthasarathi, R.; Elango, M.; Subramanian, V.; Krishnamoorthy, B. S.; Gutierrez-Oliva, S.; Toro-Labbe´, A.; Roy, D. R.; Chattaraj, P. K. J. Phys. Chem. A 2007, 111, 9130. (58) Morell, C.; Grand, A.; Toro-Labbe´, A. J. Phys. Chem. A 2005, 109, 205. (59) Morell, C.; Grand, A.; Toro-Labbe´, A. Chem. Phys. Lett. 2006, 425, 342. (60) Geerlings, P.; De Proft, F. Phys. Chem. Chem. Phys. 2008, 10, 3028. (61) Ayers, P. W.; Morell, C.; De Proft, F.; Geerlings, P. Chem.sEur. J. 2007, 13, 8240. (62) Ayers, P. W.; Parr, R. G. J. Chem. Phys. 2008, 129, 054111. (63) Breneman, C. M.; Wiberg, K. B. J. Comput. Chem. 1990, 11, 361. (64) Labet, V.; Grand, A.; Cadet, J.; Eriksson, L. A. ChemPhysChem 2008, 9, 1195. (65) Eyring, H. Chem. ReV. 1935, 17, 65.
JP9049044