Modeling the Dissociative Hydrolysis of the Natural DNA Nucleosides

Dec 29, 2009 - Two-dimensional PCM-B3LYP/6-31+G(d) potential energy surfaces for the hydrolysis of the four natural 2′-deoxyribonucleosides ...
0 downloads 0 Views 3MB Size
Download PDF
1104

J. Phys. Chem. B 2010, 114, 1104–1113

Modeling the Dissociative Hydrolysis of the Natural DNA Nucleosides Jennifer L. Przybylski and Stacey D. Wetmore* Department of Chemistry and Biochemistry, UniVersity of Lethbridge, 4401 UniVersity DriVe, Lethbridge, Alberta T1K 3M4, Canada ReceiVed: October 14, 2009; ReVised Manuscript ReceiVed: NoVember 30, 2009

Two-dimensional PCM-B3LYP/6-31+G(d) potential energy surfaces for the hydrolysis of the four natural 2′-deoxyribonucleosides (2′-deoxyadenosine, 2′-deoxyguanosine, 2′-deoxycytidine, and thymidine) are characterized using a model that includes both implicit (bulk) solvent effects and (three or four) explicit water molecules in the optimization routine. For the first time, the experimentally predicted dissociative (SN1) mechanism is found to be favored over the synchronous (SN2) pathway for all nucleosides studied. Due to the success of our model in stabilizing the charge-separated intermediates along the SN1 pathway, it is proposed that the new model presented here is the smallest system capable of generating the experimentally predicted oxacarbenium cation intermediate. We therefore stress that dissociative mechanisms should be studied with methodologies that account for the (bulk) environment in the optimization routine, where these effects are often only included as a correction to the energy in the current literature. In addition to accounting for charge stabilization through implicit solvation, nucleophile activation and leaving group stabilization should also be explicitly introduced into the model to further stabilize the system. Our work also emphasizes the importance of studying the Gibbs surface, which in some cases provides a better description of chemically important regions of the reaction surface or changes the calculated trend in the magnitude of dissociative barriers. In addition, it is proposed that the methodology presented in this study can be used to calculate uncatalyzed deglycosylation barriers for a range of DNA nucleosides, which when compared to the corresponding enzymecatalyzed reactions, will allow the prediction of the rate enhancement (barrier reduction) due to the enzyme. Introduction Nucleic acids are involved in a wide variety of physiological reactions, including depurination or depyrimidination,1,2 deamination of the nucleobase,3,4 and phosphodiester cleavage of the backbone.5 These reactions can occur spontaneously1,2,4 or be caused by metabolic factors such as radicals.5 They may also be facilitated by an enzyme. For example, glycosylases6,7 (e.g., UDG, MutY), phosphorylases8 (e.g., thymidine phosphorylase) and hydrolases9 (e.g., Ricin A-Chain) all catalyze the cleavage of the sugar-nucleobase (N-glycosidic) bond. Due to the importance of these classes of enzymes with respect to natural processes (e.g., nucleotide recycling9 and aging10) and disease (e.g., cancer11 and HIV12), their function has been the focus of many computational groups.13-34 These theoretical studies have ranged in complexity from considerations of relative substrate stability13-20 to small model studies of the deglycosylation reaction21-29 to large model examinations of enzymatic mechanisms and rates.30-34 To properly determine the rate enhancement (or catalytic efficiency) associated with an enzyme, the calculated barriers for the enzyme-catalyzed mechanism must be compared to the related spontaneous reaction in (aqueous) solution.35 In the case of glycosylases, phosphorylases, and hydrolases, the baseline solution reaction corresponds to hydrolytic cleavage of the glycosidic bond in nucleotides or nucleosides. Therefore, it is essential to carefully study the uncatalyzed hydrolysis of nucleic acids. Unfortunately, in many cases, the uncatalyzed reaction occurs so slowly that experimental measurements of the reaction are extremely difficult. As a result, experimental techniques that enhance the rate in a predictable fashion are often utilized. For * Corresponding author. E-mail: [email protected].

example, studies on the hydrolysis of 2′-deoxyribonucleosides are generally carried out in the presence of acid,36-38 or under high temperatures,2,37,39 to accelerate the reaction. Standard methodologies are available for extrapolating high-temperature results to physiological conditions with a fair degree of accuracy, and this approach has been utilized to study a broad variety of reactions involving biomolecules (for example, hydrolysis of peptides).40 However, in cases where this extrapolation cannot be performed (for example, the Arrhenius plot is nonlinear) and experiments cannot be easily carried out under physiological conditions, computational techniques can contribute to determining reaction barriers. An assortment of experimental work on the nonenzymatic deglycosylation of nucleic acids has been reported in the literature.1,2,36-39 The hydrolytic cleavage of the glycosidic bond in both nucleosides and nucleotides has been proposed to be dissociative (SN1) in nature rather than synchronous (SN2), which is supported by positive entropies of activation for these reactions.37-39 In aqueous solution, the rate-limiting step for the glycosidic bond cleavage is proposed to be dissociation of the nucleobase from the sugar to yield an oxacarbenium cation and a nucleobase anion (Figure 1). A water molecule then adds to the sugar moiety of the charge-separated intermediate, and loses a proton to the surrounding environment. Similarly, the nucleobase anion accepts a proton from the solvent. Despite experimental predictions, some computational studies have focused solely on the synchronous mechanism for nucleoside hydrolysis.13,16,27 Other studies that examine both the SN1 and SN2 mechanisms report prohibitively large activation energies for the SN1 dissociation step, and thus conclude that the synchronous mechanism is favored.25,26,29 In addition, most studies do not contain the empirically proposed oxacarbenium

10.1021/jp9098717  2010 American Chemical Society Published on Web 12/29/2009

Hydrolysis of the Natural DNA Nucleosides

Figure 1. Experimentally proposed steps for (neutral) hydrolysis of DNA nucleosides. In the first step, the base dissociates from the sugar leaving an oxacarbenium cation and nucleobase anion, which is followed by irreversible addition of a water molecule. The nucleobase product may involve a variety of tautomers, where the two most stable uracil tautomers are shown as an example.

Figure 2. Structures and atomic numbering of the 2′-deoxynucleosides considered in the present work, namely 2′-deoxyadenosine (dA), 2′deoxyguanosine (dG), 2′-deoxycytidine (dC), thymidine (dT), and 2′deoxyuridine (dU).

cation intermediate.24-26,29 For example, two studies have reported a reaction intermediate where the nucleobase reattaches to the sugar through a different connectivity compared with the corresponding natural nucleoside,26,29 such as O2 of thymine glycol bound to C1′ of the sugar moiety.29 A reaction intermediate whereby the nucleobase (directly or indirectly) abstracts a proton from C2′ of the sugar has also been reported for both nucleoside hydrolysis25,26,29 and enzymatic depurination.24 However, neither of these reaction intermediates is supported by empirical evidence. Discrepancies between experimental and computational results for nucleotide or nucleoside deglycosylation may be partially attributed to the use of gas-phase optimizations, where charge-separated systems, such as dissociative transition states, have been shown to be highly unstable.40 In all of the uncatalyzed hydrolysis computational studies mentioned above,25-27,29 and many of the enzyme-catalyzed studies,21,23,24,31,32 environmental effects (if included) were treated as a correction to the energy. Indeed, this is still the most common methodology for studying reaction mechanisms using computational chemistry in the literature, although microsolvation and large-scale models with explicit solvation are increasing in frequency.41,42 Since the environment likely has a structural role in reactions involving charge-separated intermediates, our group recently examined the hydrolysis of 2′-deoxyuridine (dU, Figure 2) using models including one to three discrete water molecules and optimizations performed in both the gas and solvent phases.28 dU was selected due to the abundance of experimental and

J. Phys. Chem. B, Vol. 114, No. 2, 2010 1105 computational literature on its deglycosylation.7,13,14,26,36-39 To determine the role of solvent on the hydrolysis mechanism, twodimensional reaction surfaces, which encompass both the experimentally predicted SN1 and computationally predicted SN2 mechanisms, were generated.28 The most commonly employed model in the literature, which contains the nucleoside and one water molecule as the reaction nucleophile, was initially examined. Gas-phase optimizations characterized an SN2 transition state, but were unable to properly describe the SN1 reaction pathway. Specifically, in the absence of leaving group stabilization, the expected charge-separated SN1 reaction intermediate collapses such that O2 of uracil is bound to C1′ or abstracts a proton from C2′. When (bulk) solvent effects are included in the optimization, the experimentally predicted oxacarbenium cation intermediate is stabilized, but the associative transition state is not properly characterized since the nucleophilic water molecule is not sufficiently activated. However, when the solvent-phase reaction surface was calculated using a model with three explicit water molecules bridging O3′ of the sugar moiety and O2 of uracil, the entire SN1 reaction pathway was revealed. Thus, it was concluded that the smallest model capable of characterizing the SN1 mechanism must include leaving group stabilization and nucleophile activation in the form of both explicit hydrogen bonds and implicit solvent effects. In order to investigate the broader applicability of the model developed in our previous study, the hydrolysis of the four natural 2′-deoxyribonucleosides, namely 2′-deoxyadenosine (dA), 2′-deoxyguanosine (dG), 2′-deoxycytidine (dC), and thymidine (dT), will be examined in the present study by generating two-dimensional reaction potential energy surfaces. Investigating the hydrolysis of the natural nucleosides will not only test the proposed hybrid solvation model for nucleoside hydrolysis but will also provide a baseline for the uncatalyzed deglycosylation barrier for comparison to future computational studies of enzyme-catalyzed deglycosylation reactions. Computational Details As discussed above, our previous study on dU hydrolysis found that gas-phase optimizations typically used in the literature to study nucleoside hydrolysis are biased toward a synchronous (SN2) mechanism. Therefore, the hydrolysis of the natural 2′deoxynucleosides will be examined using a model based on our successful dU case study.28 Specifically, nucleophile activation and leaving group stabilization will be accounted for through an explicit water network between the (water) nucleophile and (nucleobase) leaving group, and (bulk) environmental effects will be accounted for using implicit solvation in the optimization routine. As described below, the hydrolysis of the nucleosides will be investigated using potential energy surface (PES) scans, where scanning the reaction surface has proven to be a useful tool for studying other hydrolysis reactions (for example, phosphodiester cleavage43,44 and those catalyzed by glycosylases32). Furthermore, the importance of analyzing nonstationary points has recently been emphasized in the literature,44 where this analysis has been shown to play a critical role in understanding biological reactions with broad, flat reaction surfaces. As per our previous work on dU,28 all geometry optimizations were carried out using the B3LYP density functional with the 6-31+G(d) basis set. This combination of basis set and method is often utilized to study biological reaction mechanisms.13,16,19,23-29,43 A larger basis set was previously tested on the dU reaction surface, and all relative energies were within 3 kJ mol-1 of those calculated with the smaller basis

1106

J. Phys. Chem. B, Vol. 114, No. 2, 2010

set.28 Implicit (bulk) solvation effects due to the water environment (ε ) 78.39) were incorporated in the optimizations using the integral-equation formalism polarizable continuum model (IEF-PCM),45 with UAHF radii (rmin ) 0.5 and ofac ) 0.8). We note that UAHF radii were utilized since they yield similar results to UAKS radii,46 and provide better hydration energies.47 Frequency calculations were carried out in (bulk) solvent at the same level of theory to confirm that any imaginary frequencies correspond to movements along the reaction pathway, which represents the intrinsic reaction coordinate (IRC).48 The Gibbs energy for each point on the surface was calculated using a standard approach.47 Specifically, the scaled (0.9806) zero-point vibrational energy, thermal correction and (unscaled) entropy term obtained from the solvent-phase frequency calculations on each geometry were added to the corresponding solvent-phase (optimization) energy. All optimization and frequency calculations were carried out using Gaussian 03 (Rev. D.02).49 In our previous study on dU,28 the reactant structure was generated by optimizing the bare nucleoside using B3LYP/631+G(d) and adding three water molecules by inspection to connect the C3′-hydroxyl to O2 of uracil. Subsequently, Monte Carlo sampling of the water positions was carried out at 300 K using the default settings in HyperChem50 with the AMBER51 molecular mechanics force field. The 50 lowest energy structures were optimized at B3LYP/6-31+G(d). After optimization, the structure with the lowest energy (which also corresponded to the largest number of Monte Carlo parent structures) was selected as the reactant geometry for the PES scans (see below) and was reoptimized with PCM-B3LYP/6-31+G(d). Since the natural pyrimidines are structurally similar to dU, the pyrimidine reactant structures were generated by modifying the dU reactant to the correct nucleobase, and this structure was optimized to a minimum using PCM-B3LYP/6-31+G(d). The initial geometry for the purine nucleosides was designed to resemble the dU reactant geometry. However, four water molecules (rather than three) connect the C3′-hydroxyl to N3 of the purine base (Figure 2) since the larger purine ring system increases the distance of the hydrogen-bond network. In fact, it was found that a purine model containing only three discrete water molecules creates strain in the system and alters the sugar puckering compared to the nucleoside optimized in fully explicit water.52 Subsequently, as done for dU, the optimum position and orientation of the water molecules were determined using Monte Carlo simulations. Due to the additional water molecule, the lowest energy structures within a 10 kJ mol-1 range (∼80 structures) were further optimized at PCM-B3LYP/6-31+G(d). In all cases, the reactant utilized to generate the PES contained C2′-endo sugar puckering, the C5′-hydroxyl group was oriented away from the nucleobase, and the nucleobase was in an anti χ conformation. Once a reactant geometry was selected, the glycosidic bond length (C1′ · · · N1 (pyrimidine) or C1′ · · · N9 (purine)) was set to 1.5 Å, the nucleophile distance (C1′ · · · Onuc) was set to 3.4 Å, and the remainder of the system was allowed to relax at the PCM-B3LYP/6-31+G(d) level. The reaction surface was generated by incrementally elongating or compressing one of the distances by 0.2 Å and relaxing the system to the new constraints. In total, the glycosidic bond length was sampled between 1.3 and 4.3 Å, while the nucleophile distance spanned 1.2-3.8 Å. The resulting potential and Gibbs energy surfaces are plotted with the glycosidic bond length along the x-axis and the nucleophilic distance on the y-axis, such that the reactant is always found in the bottom left corner and the product in the top right corner of the plot. The relative energies are given in color with the reactant occurring at 0.0 kJ mol-1 (purple) and

Przybylski and Wetmore each 10 kJ mol-1 increment represented by a single change in color. For clarity, the resolution is increased for the transition state and intermediate regions by using 2 kJ mol-1 increments, which is indicated by black vertical bars in the legends for the contour plots. Due to the use of constrained geometries, all energies are reported as an average value approximated from a region of the surface, where a representative structure from each region is provided in the figures. Results and Discussion In our previous work, a model for 2′-deoxynucleoside hydrolysis was developed using the smallest nucleobase, uracil.28 It was found that the most effective computational model incorporates a hydrogen-bonding network between explicit water molecules and the nucleobase, as well as (bulk) solvent effects, during the optimization routine. This methodology is now applied to the natural DNA nucleosides to determine if our model is applicable to other related systems. Initially, the hydrolysis of the pyrimidines (dT and dC) using a model with three explicit water molecules will be discussed. The calculated reaction surfaces can be found in Figure 3, representative structures are given in Figure 4, and typical reaction energies estimated from the respective surfaces are provided in Table 1. The results from our previous dU study are included with the other pyrimidine data for comparison. Subsequently, a discussion of the hydrolysis of the purines (dA and dG) using a model with four explicit water molecules will be provided (Figures 5 and 6, and Table 2). Pyrimidines. Thymidine. The potential energy surface for dT hydrolysis (Figure 3A) indicates that the reaction path of slowest ascent involves dissociation of the glycosidic bond to a shallow intermediate, which is followed by nearly barrierless association of the nucleophilic water molecule before reaching the product well. A very flat dissociation surface has been reported previously for dG depurination22 and dU hydrolysis,28 where it was found that calculating the Gibbs energy surfaces enhances relatively flat regions of the corresponding PES. Indeed, the Gibbs energy surface for dT hydrolysis (Figure 3B) reveals a deeper intermediate well (∼7 kJ mol-1 more stable than the dissociative transition state, see Table 1) and amplifies the region of the plot that corresponds to the associative transition state. With the exception of the glycosidic bond and nucleophile distances, the geometrical parameters are very similar for structures selected from chemically important regions of the potential energy and Gibbs energy surfaces. Therefore, the representative geometries along the reaction pathways displayed in Figure 4 were selected from the PES, which also allows for direct comparison with our dU study and other previous literature. In the reactant (T-R, Figure 4A), the water chain interacts with O2 of the nucleobase and the C3′-hydroxyl of the sugar moiety. The dissociative transition state (T-TSD, Figure 4A) occurs near a glycosidic bond length of 2.9 Å, and involves a ∼130 kJ mol-1 energy barrier (∼96 kJ mol-1 on the Gibbs surface, Table 1). The ∆Gq for dT dissociation is only slightly smaller than that calculated for dU (∼100 kJ mol-1).28 Although the direct interaction between the C3′-hydroxyl and the nucleophile in our model may not be present in a fully (explicitly) solvated system, we anticipate the effects of removing this interaction on reactant stabilization and nucleophile activation should nearly cancel, and therefore this interaction should not have a large effect on the reported barriers. The hydrogen bond between the nucleobase and the water chain decreases from 1.818 Å in the reactant complex to 1.654 Å in the transition

Hydrolysis of the Natural DNA Nucleosides

J. Phys. Chem. B, Vol. 114, No. 2, 2010 1107

Figure 3. PCM-B3LYP/6-31+G(d) potential (A) and Gibbs (B) energy surfaces for the hydrolysis of dU (top), dT (middle), and dC (bottom) (kJ mol-1; color scale provided in legend; see Computational Details; dU data from ref 28).

state, which aids stabilization of the negative charge accumulating on the thymine moiety. The sugar group is similarly stabilized by the nucleophilic water molecule, which moves closer to the sugar and reorients to direct a lone pair toward the positive charge accumulating on C1′ and O4′. These structural changes are parallel to those previously reported for dU hydrolysis,28 where the O2 carbonyl is involved in hydrogen bonds to stabilize the charge build up on the nucleobase. The intermediate (T-I, Figure 4A) is 4 kJ mol-1 more stable than the dissociative transition structure on the PES (7 kJ mol-1 on the Gibbs surface, Table 1). In T-I, the thymine anion rotates

to introduce a weak interaction between N1 and the C1′ proton. Due to differences in O2 basicities,14 this contrasts the structure of the dU intermediate, where O2 of uracil is involved in an interaction with C1′.28 Most importantly, our dT hydrolysis model does not find either the O2-bound or C2′-H abstraction reaction intermediates reported in a gas-phase study on the hydrolysis of the related thymidine glycol,16 which emphasizes the importance of including sufficient explicit and implicit solvent in the optimization routine to remove the bias toward charge neutralization. The second, associative, transition state (T-TSA, Figure 4A) involves a very small barrier (∼1 kJ mol-1

1108

J. Phys. Chem. B, Vol. 114, No. 2, 2010

Przybylski and Wetmore

Figure 4. Representative structures for pyrimidine hydrolysis determined from the PES for dU (top), dC (middle), and dT (bottom), including the SN1 pathways (A) and SN2-like transition states (B). (Select PCM-B3LYP/6-31+G(d) distances (Å) provided including glycosidic and nucleophilic distance represented as [C1′ · · · N1,C1′ · · · Onuc], dU data from ref 28.)

TABLE 1: Relative Potential and Gibbs Energies Estimated for the Dissociative and Concerted Hydrolysis of the Pyrimidine Nucleosidesa,b R

TSD

dUc dT dC

0 0 0

120 130 141

dUc dT dC

0 0 0

100 96 110

I ∆E 119 126 139 ∆G 90 89 96

TSA

P

TSC

123 127 142

30 38 27

160 150 167

115 102 122

35 28 17

140 140 148

a

All PCM-B3LYP/6-31+G(d) relative energies reported with respect to the corresponding reactant. b See Figure 3 for contour plots and Figure 4 for representative structures. c Reference 28.

on the PES; ∼13 kJ mol-1 on the Gibbs surface). The associated imaginary frequency corresponds to the synchronous motion of the nucleophilic water molecule toward C1′ and stretching of the nucleophilic water O-H bond involved in the hydrogenbonding network to the nucleobase, which is the initial stage of proton transfer to the leaving group. Indeed, in the product (T-P, Figure 4A), this proton completely transfers to the neighboring water molecule, and a proton shuffle occurs through the water network to O2 of the thymine anion to yield the O2 protonated thymine tautomer. A similar proton shuffle has also been reported for thymidine phosphorylase, where a salt bridge transfers a proton from the nucleophile to O2 of thymine.31 The total reaction energy for the hydrolysis is ∼38 kJ mol-1 (∆Gq ≈ 28 kJ mol-1). As seen previously in our study of dU hydrolysis,28 a transition structure corresponding to a synchronous mechanism can be found on the surface (T-TSC, Figure 4A) with a ∆Eq of ∼150 kJ mol-1 (∆Gq of ∼140 kJ mol-1, Table 1). Therefore, our new model correctly predicts that the SN1 mechanism is more favorable than the SN2 mechanism. The synchronous

transition structures for dT and dU (U-TSC, Figure 4A) are structurally very similar and the corresponding barriers are within 10 kJ mol-1 (Figure 3A). Our group recently studied the hydrolysis of the 2′-deoxynucleotides in the gas phase using full optimizations, a nucleotide model with methoxyl groups at C3′ and C5′ (rather than hydroxyl groups), and a formate-water complex as the nucleophile.27 This study found a (B3LYP/6311+G(2d,p)//B3LYP/6-31G(d)) hydrolysis barrier of 112.6 kJ mol-1 when (bulk) solvent effects were included as an energy correction. Since the formate-water nucleophile is more reactive than the one used in the present study, the previously calculated barrier corresponds well with the present work. Therefore, despite a large difference in the structures characterized for the SN1 pathway between previous (gas phase) work and the current study, the concerted transition state is very similar to previous (gas phase) literature. This suggests that inclusion of solvent effects during the optimization step is more important for dissociative than synchronous structures, which corresponds with previous literature findings that the environment is important for correctly predicting charge-separated structures.53 In the Introduction, it was proposed that the methodology presented in this work could be used as a baseline for enzymecatalyzed deglycosylation studies. The example of human thymidine phosphorylase, which catalyzes the depyrimidination of the thymidine nucleoside, will be utilized to illustrate this application since reliable experimental data and computational studies of the enzymatic reaction are available in the literature. In a previous computational study by Rick et al.,31 the energy of activation for the enzyme-catalyzed dT depyrimidination was calculated to be 55.2 kJ mol-1. In the present work, the activation energy for dT hydrolysis is found to be ∼130 kJ mol-1, which indicates that thymidine phosphorylase lowers the deglycosylation barrier by ∼75 kJ mol-1. When the catalyzed (kcat ) 8.2 s-1)54 and uncatalyzed (k25°C ) 2.5 × 10-10 s-1)39 experimental rates of thymidine deglycosylation are compared,

Hydrolysis of the Natural DNA Nucleosides

J. Phys. Chem. B, Vol. 114, No. 2, 2010 1109

Figure 5. PCM-B3LYP/6-31+G(d) potential (A) and Gibbs (B) energy surfaces for the hydrolysis of dA (top) and dG (bottom) (kJ mol-1; color scale provided in legend; see Computational Details).

Figure 6. Representative structures for purine hydrolysis determined from the PES for dA (top) and dG (bottom), including the SN1 pathways (A) and SN2-like transition states (B). (Select PCM-B3LYP/6-31+G(d) distances (Å) provided including glycosidic and nucleophilic distance represented as [C1′ · · · N9,C1′ · · · Onuc].)

a rate enhancement of 3.3 × 1010 is calculated, which corresponds to a decrease in the Gibbs energy of activation by ∼60 kJ mol-1. This value is very close to the change in activation energy found by comparing the calculated barriers.55 Although different methods were implemented for the catalyzed (B3PW91) and uncatalyzed (B3LYP) barriers, this example illustrates how the results presented in this work can act as a baseline for calculating enzymatic rate enhancements. In addition, the methodology described here may provide a useful approach for

estimating enzymatic rate enhancements when the experimental data for the corresponding uncatalyzed barriers is not available. 2′-Deoxycytidine. When the PES for the hydrolysis of the pyrimidines are compared (Figure 3A), the most striking difference is that most of the PES for dC (orange) is elevated in energy in comparison to dU and dT (green) relative to the corresponding reactants. Indeed, the dissociative ∆Eq estimated for dC is ∼141 kJ mol-1, which is 10-21 kJ mol-1 higher than the corresponding barrier for dT and dU. The large dissociation

1110

J. Phys. Chem. B, Vol. 114, No. 2, 2010

Przybylski and Wetmore

TABLE 2: Relative Potential and Gibbs Energies Estimated for the Dissociative and Concerted Hydrolysis of the Purine Nucleosidesa,b R

TSD

I

TSA

P

TSC

dA dG

0 0

134 138

∆E 128 132

133 136

0 38

144 150

dA dG

0 0

114 108

∆G 96 98

120 120

11 28

140 158

a All PCM-B3LYP/6-31+G(d) relative energies reported with respect to the corresponding reactant. b See Figure 5 for contour plots and Figure 6 for representative structures.

barrier can be attributed to lower stability of the cytosine anion compared to the uracil and thymine anions, which mimics the relative N1 acidities of the conjugate acids.14,56 The corresponding Gibbs energy surface for dC hydrolysis is presented in Figure 3B (bottom), where the ∆Gq for dissociation is also larger than dU and dT by 10-14 kJ mol-1 (Table 1). This agrees with experimental data where the Gibbs energy of activation for dC hydrolysis is larger than for hydrolysis of dU or dT, where both the measured and calculated differences in the pyrimidine reaction barriers are small.39 Similar to the dissociative transition state reported for dU and dT, the transition state for cytidine deglycosylation (C-TSD, Figure 4A) contains a strong hydrogen bond between O2 of the nucleobase and the neighboring water molecule (1.679 Å). From the dissociative region, the ∆Gq for the SN1 mechanism is estimated to be 110 kJ mol-1. In the intermediate (C-I, Figure 4A), the cytosine anion is additionally stabilized by a C1′sH · · · O2 hydrogen bond (1.883 Å), which is a tighter contact than seen for either dU28 or dT. This may arise since cytosine is more stable when protonated at O2 than N1 (the canonical form).57 Therefore, the cytosine anion interacts stronger via O2 than N1, and the C1′sH · · · N1 interaction seen for dT is not present for dC. After a 3 kJ mol-1 associative barrier (26 kJ mol-1 on the Gibbs surface), the product (C-P, Figure 4A) is found, which yields a reaction energy of ∼27 kJ mol-1 (∆G of ∼17 kJ mol-1). Both the associative transition structure and product are structurally similar to the intermediate, involving transfer of a proton from the nucleophilic water molecule to O2 of the cytosine anion with very little rotation of the nucleobase. A synchronous pathway can also be followed on the dC hydrolysis reaction surface through a concerted transition state (C-TSC, Figure 4B), which has an associated ∆Eq of ∼167 kJ mol-1 (∆Gq of ∼148 kJ mol-1). Analogous to dU and dT, the energy of activation for the concerted (SN2) pathway is larger than for the dissociative (SN1) pathway on both the PES and Gibbs surface, which indicates that our computational model predicts the correct mechanism regardless of the surface analyzed. In addition, the calculated concerted transition structure resembles those found for the other pyrimidines, where a strong, discrete hydrogen bond to O2 stabilizes the nucleobase leaving group. This is similar to the SN2 study by our group,27 where the fully optimized transition states are very similar for all pyrimidines examined. Summary. For the first time, the experimentally predicted dissociative (SN1) mechanism37-39 has been calculated to be lower in energy than the synchronous (SN2) pathway for the hydrolysis of the N-glycosidic bond in the pyrimidine DNA nucleosides. In addition, the general trend in the calculated Gibbs energies of activation for pyrimidine hydrolysis (dT < dU ,

dC, Table 1) agrees with the empirical finding that dC hydrolysis is slower than the hydrolysis of either dU or dT, where both the measured and calculated difference in ∆Gq across the pyrimidines are small.39 Furthermore, the rate enhancement derived experimentally for thymidine phosphorylase can be reasonably reproduced through comparison of our uncatalyzed barrier with a previously calculated enzyme-catalyzed barrier, which indicates that rate enhancements for less well-studied systems can be accurately predicted using calculations. The present work is, to the best of our knowledge, the first time that these reactions have been examined using a model that includes (bulk) solvent effects in the optimization routine and we attribute the successes of our model to this treatment of the environment. Indeed, when solvent-phase optimizations are performed, the reaction intermediates previously seen in the gas phase that involve the nucleobase anion either bonding to, or abstracting a proton from, the sugar moiety are not obtained. Since the corresponding high-energy reaction intermediates have also been previously reported in computational studies of purine deglycosylation,24 solvent-phase reaction scans will subsequently be implemented to describe dA and dG hydrolysis and thereby determine if these reaction mechanisms can be accurately predicted using the same methodology. Purines. 2′-Deoxyadenosine. As mentioned in the Computational Details, the dA hydrolysis model includes an additional water molecule in the hydrogen-bonding network between the sugar moiety and the nucleobase compared with the pyrimidine models, which prohibits an accurate comparison of the calculated deglycosylation barriers for the purines and pyrimidines. The potential and Gibbs energy surfaces generated using our dA hydrolysis model are given in Figure 5, A and B (top). As found for the pyrimidines and predicted experimentally,39 the path of slowest ascent follows a dissociative mechanism on both surfaces (Table 2). The PES for dA has a better defined (SN1) intermediate region compared with the pyrimidines (Figure 3A). A small associative barrier is seen before the product is formed, where this barrier is larger on the Gibbs surface than on the potential energy surface. In the reactant (A-R, Figure 6A), the adenine amino group is planar, and the water hydrogen-bonding network interacts with N3 of the nucleobase. The dissociative transition state (ATSD, Figure 6A) occurs at a ∆Eq of ∼134 kJ mol-1 (∆Gq ≈ 114 kJ mol-1) (Table 2). The dissociation is spontaneous and does not involve protonation of the nucleobase as discussed in other studies of dA deglycosylation.13,19 In the intermediate (AI, Figure 6A), the amino group begins to pucker, where the protons are directed toward the O4′ side of the nucleobase. As seen for the pyrimidines, the water-nucleobase hydrogen bond tightens at extended glycosidic bond lengths with the average distance being ∼1.78 Å for adenine. In addition, there is a weak interaction between the C1′ proton of the sugar moiety and N9 of the adenine anion, which is similar to the intermediate found for dT (T-I, Figure 4A). The intermediate exists in a < 5 kJ mol-1 well on the PES (