Designing an Appropriate Computational Model for DNA Nucleoside

Apr 9, 2009 - In the case of 2′-deoxyuridine (Figure 1), the damaged base is most ... dU is an excellent case study due to the abundance of experime...
0 downloads 0 Views 4MB Size
Download PDF
J. Phys. Chem. B 2009, 113, 6533–6542

6533

Designing an Appropriate Computational Model for DNA Nucleoside Hydrolysis: A Case Study of 2′-Deoxyuridine Jennifer L. Przybylski and Stacey D. Wetmore* Department of Chemistry and Biochemistry, UniVersity of Lethbridge, 4401 UniVersity DriVe, Lethbridge Alberta T1K 3M4 Canada ReceiVed: NoVember 28, 2008; ReVised Manuscript ReceiVed: March 10, 2009

This study uses a variety of computational models and detailed, systematic potential energy surface scans to examine the hydrolysis of 2′-deoxyuridine. First, the unimolecular cleavage was studied using a model that only includes the nucleoside. Although comparison of experimental and (PCM-B3LYP/6-31+G(d)) calculated (Gibbs energy) barriers confirms that hydrolysis occurs via a fully dissociative (SN1) mechanism with a ratelimiting step of glycosidic bond dissociation, this model does not provide a complete picture of the hydrolysis mechanism. When the model is expanded to include one explicit water nucleophile, gas-phase optimizations are unable to model charge separation in the reaction intermediate, while optimizations that implicitly incorporate the effects of bulk solvent do not accurately model the second reaction step (nucleophilic attack following dissociation) due to insufficient (water) nucleophile activation and (uracil anion) leaving group stabilization. Further expansion of the model to include three explicit water molecules allows for discrete proton transfer from the water nucleophile to the uracil anion, and thereby generates smooth reaction surfaces for both the dissociative (SN1) and concerted (SN2) pathways. Furthermore, for the first time, this computational model for the uncatalyzed hydrolysis of the N-glycosidic bond in a nucleoside predicts that the dissociative mechanism is more favorable than the concerted pathway, which supports experimental findings. It is also found that although (implicit) solvent-phase single-point calculations on gas-phase geometries can yield similar energies to solvent-phase optimizations, the geometries can be very different and not all potential reaction routes can be fully characterized. Therefore, care must be taken when interpreting mechanistic information obtained from gas-phase structures. This work provides a template for generating other nucleoside or nucleotide hydrolysis models including those relevant to both uncatalyzed and enzyme-catalyzed reactions. Introduction Deoxyribonucleic acid (DNA) uses sequences of four 2′deoxyribose nucleotides to store all information required by a cell to survive and proliferate. Although DNA is very stable, it is susceptible to internal and external attack, which initially leads to relatively small nucleobase modifications and subsequently causes larger problems such as cell mutation or death.1-5 For example, the natural base cytosine is deaminated to uracil (an RNA base) upward of 200 times per cell per day, which results in stable GC to AT point mutations after replication6,7 and has been linked to cancer.8 Therefore, nature has developed a series of repair mechanisms to ensure genome stability. In the case of 2′-deoxyuridine (Figure 1), the damaged base is most commonly eliminated through the base excision repair (BER) pathway,9 where uracil DNA glycosylase (UDG) uses a water molecule as a nucleophile to remove uracil from the sugar-phosphate backbone6,7,10-12 and other enzymes complete the repair process by replacing the base-free sugar with the correct nucleotide.6,9 Since UDG was the first glycosylase discovered,10 it has been the focus of many experimental research groups. However, relatively few computational analyses have considered its mechanism of action. The computational work published to date can be classified into three categories: (1) small model transition state studies that recreate experimentally measured kinetic isotope effects,13 (2) small model studies on the glycosidic bond * To whom correspondence [email protected].

should

be

addressed.

E-mail:

Figure 1. Structure and common atomic numbering for 2′-deoxyuridine (dU).

stability of 2′-deoxyuridine,14,15 and (3) large model hybrid quantum mechanics/molecular mechanics (QM/MM) mechanistic studies.16 Each computational analysis has focused on a different aspect of the reaction catalyzed by UDG and has employed a variety of approaches to reduce or divide the complexity of the computational model depending on the overall goal of the work. As a result, it is difficult to compare the results from these accounts. Furthermore, since each group came to diverse (and sometimes opposing) conclusions, many mechanistic unknowns surrounding UDG still exist. We believe that before this enzymatic system can be properly modeled and thereby understood, an accurate computational model that reproduces the experimentally predicted reaction mechanism for the nonenzymatic (uncatalyzed) hydrolysis reaction must be identified. Because of the important cellular role of hydrolases, which are responsible for ribonucleotide recycling, the majority of computational work on N-glycosidic bond cleavage in the literature has focused on ribonucleotide hydrolysis.17,18 However,

10.1021/jp810472q CCC: $40.75  2009 American Chemical Society Published on Web 04/09/2009

6534

J. Phys. Chem. B, Vol. 113, No. 18, 2009

it is difficult to gain information about deoxyribonucleotide hydrolysis using ribose models since the additional (C2′) hydroxyl group changes the hyperconjugation and steric roles of the sugar moiety, and therefore both the mechanism and the kinetics for glycosidic bond cleavage may be altered.19 Thus far, very few computational studies have examined the glycosidic bond cleavage in 2′-deoxyribose nucleosides,15,20-24 and most published work has focused on a particular enzymatic mechanism.20-22 Nevertheless, in some cases the computational models are small enough to be considered nonenzymatic systems. For example, many of the reaction schemes only include the nucleobase, sugar moiety and nucleophile, and are devoid of specific active-site amino acid residues.21b,22 Furthermore, previous computational surveys of DNA deglycosylation reactions have considered a variety of nucleophiles (e.g., phosphate20 or amino21 groups), where studies highlighting nucleotide hydrolysis are sparse in the literature.22,23 Rather than varying the nucleophile, other studies have modified the nucleobase through, for example, protonation (acid-catalyzed hydrolysis)23 or platination.24 In addition, most of the hydrolysis literature focuses on depurination20-24 of 2′-deoxyribonucleotides rather than depyrimidination,15 and experiments suggest that hydrolysis of the purines and pyrimidines occur at rates differing by orders of magnitude.25-27 The above discussion reveals that the recent computational literature on DNA deglycosylation has left two important areas unaddressed: (1) cleavage within pyrimidine nucleosides or nucleotides and (2) the uncatalyzed (nonenzymatic) hydrolysis mechanism. Consequently, the Wetmore group previously investigated the hydrolysis of the smallest pyrimidine nucleoside, 2′-deoxyuridine (dU), which, as discussed above, is a wellstudied form of DNA damage. A variety of mechanisms were considered, including unimolecular cleavage and hydrolysis by water or (partially) activated water.15 While this study gave important insight into the effect of nucleophile strength on the reaction kinetics, as well as the role of interactions between the nucleobase and the surrounding environment, the resulting geometries and energies are not compatible with experimental data25-27 for the uncatalyzed hydrolysis. Similar to other computational studies of the hydrolysis of nucleosides and nucleotides,20-24 this previous computational work15 used very small model systems containing only one water molecule as the nucleophile, which was partially activated through hydrogen bonds to another (model amino acid) molecule, and did not account for (bulk) solvent effects. However, (bulk) solvent effects may be important and previous literature on small organic molecules indicates that some reactions (such as proton transfer,28 tautomerization,29 and deamination30) require additional explicit water molecules.31 Therefore, further work on the nonenzymatic hydrolysis of 2′-deoxyuridine is required to understand how to accurately model this reaction. The present study focuses on systematically identifying an appropriate computational model for the hydrolysis of the 2′deoxyuridine nucleoside. dU is an excellent case study due to the abundance of experimental work on its hydrolysis,25,26 including recent Gibbs energy data for the nucleoside,27 which will permit accurate benchmarking of our computational results. We use a combination of cluster (explicit solvent molecules) and continuum (implicit solvation) techniques to vary the degree of solvation and thereby systematically increase the complexity of the computational model. The cluster-continuum approach has been previously used to study solvent effects32 and has been shown to provide a fair approximation for protic solvent effects on reaction mechanisms. We will begin our study by considering

Przybylski and Wetmore the most simplistic computational model that includes one explicit solvent molecule since this model has been shown to be unacceptable for other systems,30,31 but is commonly used in the literature to study DNA hydrolysis reactions.15,22,23 The effects of the bulk (continuum) solvent will be considered in single-point calculations (to determine the energetic effects of the solvent) and optimizations (to determine the structural effects). This varied approach is important since, although solvation effects are routinely incorporated post-optimization in studies of enzymatic reactions,33 inclusion of (implicit) solvation during optimizations can alter the geometries of reaction intermediates due to charge stabilization.34 Through our approach, key concepts about how the 2′-deoxyuridine hydrolysis reaction can be accurately modeled will be revealed, which will have important implications for future computational research on nonenzymatic hydrolysis of other 2′-deoxyribose nucleosides as well as nucleoside hydrolysis reactions catalyzed by DNA repair and other enzymes. Computational Details Although experiments suggest that the uncatalyzed hydrolysis of the N-glycosidic bond in 2′-deoxynucleosides or nucleotides occurs via a dissociative (SN1) pathway,25-27 computational studies of these reaction mechanisms generally model hydrolysis as a concerted (SN2) reaction and do not discuss the possible SN1 pathway.22 Indeed, one study explicitly addresses the theoretical difficulties involved in modeling the dissociation step.24 To ensure all possible reaction routes are considered, the hydrolysis of 2′-deoxyuridine (dU) is systematically examined in the present work using fixed coordinate optimizations to generate detailed reaction potential energy surfaces (PES). All geometry optimizations were carried out using the B3LYP density functional method with the 6-31+G(d) basis set. Test calculations with a larger (6-31+G(d,p)) basis set on representative structures indicate that our reported relative energies change by less than 1.5 kJ mol-1 with the exception of structures with multiple proton transfers, where the relative energies change by less than 3 kJ mol-1. Indeed, previous work in the Wetmore laboratory also tested a variety of basis sets (including 6-31G(d) and 6-31+G(d,p)) and found no significant differences in the geometries or barriers for nucleoside hydrolysis.15 The method and basis set implemented in the present work has also been recently used in potential energy surface scans for hydrolysis of a phosphodiester,35 and B3LYP has also been widely used to study nucleic acid glycosidic bond cleavage,15,20b,21a,b,22,23 which will allow direct comparison to our results. Prior to performing the PES scan, a variety of conformations of the dU nucleoside were optimized to identify the lowest energy structure. Subsequently, discrete water molecules were introduced by inspection and Monte Carlo calculations were conducted to sample possible water positions. The most highly populated structure following B3LYP/6-31+G(d) optimizations was used to study the PES for hydrolysis of dU. In all cases, the glycosidic bond length (N1-C1′) in the reactant complex was ∼1.5 Å and the nucleophile distance (Ow-C1′) was ∼3.4 Å. To scan the PES, the glycosidic bond length was varied between 1.3 and 4.3 Å, while the nucleophile distance was varied between 1.2 and 3.8 Å. Both distances were held fixed at 0.2 Å increments and therefore step-by-step hydrolysis surfaces containing ∼200 individually optimized structures were obtained. Implicit (bulk) solvation effects due to water were incorporated using the integral-equation formalism polarizable continuum model (IEF-PCM).36 Two separate methods for incor-

A Case Study of 2′-Deoxyuridine

J. Phys. Chem. B, Vol. 113, No. 18, 2009 6535

porating bulk solvation effects were utilized: (1) solvent-phase (PCM) single-point (energy) calculations on gas-phase geometries and (2) (implicit) solvent-phase (PCM) optimizations. In both cases, the same level of theory employed in the gas-phase optimizations was implemented. A correction for potential basis set superposition error was not included in the relative energies since difficulties defining the counterpoise fragments for nonminimum structures (i.e., transition states) are anticipated to result in greater errors in our reported relative energies.37 The Gibbs energy surfaces were calculated in the standard fashion.38 Specifically, the surfaces generated from the gas-phase optimized geometries were calculated by adding the scaled (0.9806) zero-point vibrational energy, thermal correction and (unscaled) entropy term obtained from the gas-phase frequency calculation to the solvent-phase single-point energy. Similarly, the surfaces generated from the solvent-phase optimized geometries were calculated by adding the (scaled) zero-point energy, thermal correction, and (unscaled) entropy term calculated in water to the solvent-phase (optimization) energy. All calculations were carried out using Gaussian 03.39 Comprehensive computational details are provided in the Supporting Information. Results and Discussion High-level computational studies of enzymatic reactions typically involve reduction to a small model system.33 However, before reducing the enzymatic problem, it is important to understand how to design an accurate nonenzymatic model since features required for the nonenzymatic reaction should also be compulsory for the catalytic model. This question is addressed in the present work by considering the specific case of the uncatalyzed hydrolysis of 2′-deoxyuridine, which will allow accurate future studies of the hydrolysis of other DNA nucleosides and nucleotides, as well as hydrolysis reactions catalyzed by DNA repair (or other) enzymes. Initially, a model that only includes the nucleoside will be considered. Subsequently, models that incorporate (one to three) discrete water molecules will be investigated, where the results for each system of increasing size are discussed in detail in the subsequent sections. UniMolecular Glycosidic Bond Cleavage in dU. Experimental studies suggest that the hydrolysis of dU proceeds via a dissociative SN1 pathway with a rate-limiting step of glycosidic bond cleavage to produce a fully formed uracil anion and oxacarbenium cation.25-27 This indicates that the simplest computational model that may reproduce the experimental barrier for dU hydrolysis will correspond to the first, dissociative step of the SN1 mechanism. Thus, we initially modeled the glycosidic bond cleavage in dU in the absence of the (water) nucleophile (i.e., unimolecular cleavage). Figure 2A plots the relative energy (∆E) of the nucleoside model as a function of the glycosidic bond length and reveals that the gas-phase surface is unstable. Specifically, once the glycosidic bond is elongated to approximately 3 Å, a sharp energy reduction occurs from roughly 160 to 90 kJ mol-1. This energy drop is associated with a drastic geometrical change, where the uracil anion abstracts a proton from C2′ of the sugar moiety (Figure 2A, inset). Despite the lack of experimental evidence for a C2′-proton (C2′H) abstraction product, previous computational studies of glycosidic bond cleavage reactions in dU,15 as well as other DNA nucleosides,21b,23 which implement full (unconstrained) transition state optimizations, report similar findings. Indeed, our dU dissociation barrier estimated from the energy plateau preceding the geometrical change (∼160 kJ mol-1) is comparable to the previously published barrier (155 kJ mol-1) based on full optimizations.15,40

Figure 2. B3LYP/6-31+G(d) potential energy surface for dU unimolecular glycosidic bond cleavage calculated in the gas phase (A, blue diamonds), in the solvent phase (water) using gas-phase geometries (A, green squares), and in the solvent phase (water) using solventphase geometries (B). Representative structures along the reaction pathway are inset.

Because of structural discrepancies between the calculated gas-phase unimolecular pathway and experimental data, it is clear that gas-phase optimizations are inadequate for this system even though they are often used in the literature. Indeed, the C2′H-abstraction discussed above likely arises due to the inability of gas-phase calculations to model the charge separation present in the expected product. This shortcoming of gas-phase calculations may also contribute to the previous incorrect computational prediction that dU hydrolysis occurs via an SN2, rather than an SN1, pathway.15 Therefore, the next step to improve the computational model is to accommodate charge separation by including bulk solvent effects. Figure 2 compares the potential energy surfaces for dU deglycosylation calculated upon inclusion of bulk solvent effects through single-point calculations, which determine the energetic effects, and geometry optimizations, which determine the structural effects. Although solvent-phase single-point calculations stabilize the gas-phase activation energy by about 15 kJ mol-1 (Figure 2A), problems with the gas-phase geometry (i.e., C2′H-abstraction in the product complex) can clearly not be rectified. In contrast, solvent-phase optimizations do not lead to C2′H-abstraction in the product complex. Instead, the reaction proceeds to a charge-separated product, where N1 of the uracil anion hydrogen bonds with C1′ of the cationic sugar moiety (Figure 2B, inset). Thus, even though the barrier (∼155 kJ mol-1) for the solvent-phase optimized reaction (estimated from

6536

J. Phys. Chem. B, Vol. 113, No. 18, 2009

Przybylski and Wetmore

Figure 3. PCM-B3LYP/6-31+G(d) Gibbs energy surface for dU unimolecular glycosidic bond cleavage calculated in (bulk) water using gas-phase geometries (green diamonds) and solvent-phase geometries (red circles).

the energy plateau) is similar to the gas-phase estimated barrier (∼160 kJ mol-1), it is imperative to consider solvent-phase optimizations to obtain realistic structures along the reaction coordinate. It is also noted that the continuous increase in the solvent-phase potential energy as the glycosidic bond is elongated (Figure 2B) is likely the reason previous studies have highlighted difficulties modeling nucleoside (or nucleotide) deglycosylation reactions using full optimizations.24 Since the evaluation of Gibbs energies should permit a direct comparison of calculated and experimental barriers, the Gibbs energy profile for dU deglycosylation in water was determined using both the gas and solvent-phase optimized geometries (Figure 3). The gas-phase structures yield a flat Gibbs energy plateau at approximately 120 kJ mol-1, while the true solventphase Gibbs energy surface contains a similarly broad maximum near 115 kJ mol-1. Interestingly, considering the size of the computational model, which includes the nucleoside in the absence of a discrete (water) nucleophile, and despite difficulties associated with the accurate calculation of Gibbs energies, which are underestimated due to over estimations of the entropy term,41 the calculated ∆G‡ values are close to the most recently reported experimental value (127 kJ mol-1).27 Thus, our calculations support experimental suggestions that the hydrolysis of 2′deoxyuridine follows a dissociative, SN1, mechanism.25-27 However, while both solvation methods yield similar energies, the gas-phase optimizations breakdown at extended glycosidic bond lengths. This suggests that previous work that includes solvent effects by performing single-point calculations on gasphase geometries is likely not accurately modeling nucleoside deglycosylation reactions. dU Hydrolysis Modeled with One Discrete Water Molecule. Although studying the unimolecular glycosidic bond cleavage in dU gives valuable information about the first, dissociative, step of the hydrolysis reaction, this model does not afford a complete picture of the reaction mechanism. Therefore, the next step in systematically increasing our computational model is to include a discrete water molecule as the reaction nucleophile. While the previous section revealed the importance of including the bulk solvent in optimizations of reaction intermediates for the unimolecular bond cleavage, it is not clear how the presence of a discrete nucleophile will affect this conclusion. Furthermore, one of the goals of the present study is to determine both the energetic and structural

Figure 4. B3LYP/6-31+G(d) potential energy surface (kJ mol-1, color scale provided in legend) for the hydrolysis of dU with one explicit water molecule in the gas phase (A), in the solvent phase (water) using gas-phase geometries (B), and in the solvent phase (water) using solvent-phase geometries (C).

effects of bulk solvent. Hence, the potential energy surface for hydrolysis using a model that includes one discrete water nucleophile was initially characterized in the gas phase, and subsequently bulk solvent effects were incorporated implicitly through either single-point calculations or optimizations. Figure 4A displays the three-dimensional potential energy surface for the gas-phase hydrolysis of dU as a contour plot where the x-axis represents the distance between N1 of uracil and C1′ of the sugar (i.e., the glycosidic bond length), the y-axis represents the distance between the oxygen in the nucleophile and C1′ of the sugar (i.e., the nucleophile distance), and the relative energy is plotted in color. The top-left corner of the surface with glycosidic bond lengths between 1.3 and 2.1 Å

A Case Study of 2′-Deoxyuridine

Figure 5. Representative structures for the dissociative, SN1, (A) and concerted, SN2, (B) pathways from the gas-phase potential energy surface (Figure 4A) for dU hydrolysis modeled with one discrete water molecule. (Select B3LYP/6-31+G(d) distances (Å) provided including glycosidic bond length and nucleophile distance represented as [C1′ · · · N1,C1′ · · · Ow].)

and nucleophile distances between 1.2 and 1.8 Å was not modeled since these geometrical constraints yield highly compressed structures, and preliminary calculations (not shown) indicate that this region is much higher in energy than the rest of the surface. Representative structures throughout the reaction coordinate are displayed in the figures and were selected based on relative energy and the number of negative (imaginary) frequencies. Although structures reported as transition states are not always true first-order saddle points, the largest (and generally only) imaginary frequency corresponds to the desired motion. The reactant is the lowest energy point on the potential energy surface (purple region), where the water nucleophile hydrogen bonds to both the C3′-hydroxyl of the sugar and O2 of uracil (R, Figure 5). A dissociative SN1 transition state (denoted as TSD in Figure 5A) corresponding to uracil departing from the sugar moiety can be found, where the water nucleophile follows the uracil anion away from the sugar due to a strong Ow-H · · · O2 hydrogen bond. This hydrogen bond prevents the uracil moiety, as well as the lone pairs on the nucleophile, from efficiently stabilizing the positive charge developing on the sugar. Consequently, as the glycosidic bond is further lengthened on the gas-phase PES, unrealistic strains on the system are introduced that alleviate themselves through spontaneous geometrical rearrangements, which give rise to significantly different structures and energies (bottom right quadrant). Specifically, the red region (glycosidic bond lengths ranging between 2.9 and 3.9 Å) corresponds to the previously discussed C2′Habstraction dissociation product, while the blue energy well (glycosidic bond lengths greater than 2.9 Å and nucleophile distances ranging between 2.6 and 3.8 Å) corresponds to a reaction intermediate where O2 of uracil is fully bound to C1′ of the sugar moiety (O2-Bound, Figure 5A). Interestingly, the discontinuous region occurs at a glycosidic bond length (>2.9 Å) comparable to that discussed for the unimolecular gas-phase reaction (Figure 2A), and similarly may arise due to the instability of charge-separated species in the gas phase. Despite the collapse of the gas-phase potential energy surface, our calculations suggest that the dissociation barrier will be greater

J. Phys. Chem. B, Vol. 113, No. 18, 2009 6537 than 160 kJ mol-1, which is very similar to the value calculated for unimolecular cleavage. Since both SN1 intermediates discussed above result in a collapse of the potential energy surface, neither structure (O2bound nor C2′H-abstraction) directly leads to a hydrolysis product. However, a continuous concerted pathway (Figure 5B) can be followed from the reactant to the product, which occurs in the top right purple well on the PES. In the product complex, the uracil anion abstracts a proton from the original (water) nucleophile to form the O2 protonated uracil tautomer, which hydrogen bonds to the new C1′-hydroxyl group using both N1 and the new O2 proton (P, Figure 5B). The SN2 transition state (TSC, Figure 5B) leads to an approximate barrier of 175 kJ mol-1. Therefore, even though only the SN2 pathway can be properly characterized on the gas-phase optimized surface, the path of slowest ascent corresponds to SN1 dissociation. Since the above discussion reveals that the major problems with the gas-phase hydrolysis surface remain even when a discrete (water) nucleophile is included in the computational model, bulk solvent effects were incorporated as discussed for the unimolecular mechanism. First, solvent-phase single-point calculations were performed on gas-phase geometries. These calculations reveal that the energetic effect of the bulk solvent on the dissociative energy barrier (Figure 4B) is over 30 to 125 kJ mol-1, which is much greater stabilization than seen for the unimolecular mechanism (approximately 15 kJ mol-1, Figure 2A). When the corresponding Gibbs energy surface is calculated (Figure 6A), the dissociative ∆G‡ is approximately 100 kJ mol-1, which is roughly 20 kJ mol-1 lower than that calculated for unimolecular cleavage (120 kJ mol-1, Figure 3) and almost 30 kJ mol-1 lower than the experimental value (127 kJ mol-1).27 Most importantly, it must be emphasized that while solventphase single-point calculations on gas-phase geometries are able to provide energies that are more experimentally relevant, they do not address difficulties with the structures of reaction intermediates (i.e., the C2′H-abstraction and O2-bound states). Thus, although many groups justify gas-phase reaction mechanisms by comparing experimental barriers to calculated barriers adjusted using solvent-phase single-point calculations, this study shows that caution should be taken when interpreting mechanistic conclusions drawn from gas-phase geometries, especially for reaction pathways involving charge separation. Second, solvent-phase optimizations were performed. The resulting potential energy surface (Figure 4C) is quite different from the surfaces generated using gas-phase geometries (Figures 4A and 4B), which suggests the solvent has a large structural effect on this reaction. On the solvent-phase surface, the reactant is still the lowest energy point (R, Figure 7), and resembles the gas-phase optimized structure (R, Figure 5). However, the path of slowest ascent now leads to a dissociative intermediate, where the uracil anion is involved in a C1′-H · · · O2 hydrogen bond (I, Figure 7A). The corresponding dissociative barrier is approximately 130 kJ mol-1, where the transition structure (TSD, Figure 7A) shows partial dissociation of the glycosidic bond and formation of the new (C1′-H · · · O2) hydrogen bond observed in the intermediate. In contrast to the gas-phase optimized geometry (Figure 5A), the water nucleophile stays below C1′ of the sugar in the transition state, and thus helps stabilize the oxacarbenium cation being formed. A high energy associative transition structure can also be found (TSA, Figure 7A) with a barrier of 5 kJ mol-1 relative to the intermediate (138 kJ mol-1 relative to the reactant). However, the path of slowest ascent (dark green band, Figure 4C) does not point toward the associative transition state, but rather suggests further

6538

J. Phys. Chem. B, Vol. 113, No. 18, 2009

Figure 6. PCM-B3LYP/6-31+G(d) Gibbs energy surface (kJ mol-1, color scale provided in legend) for dU hydrolysis modeled in bulk water with one explicit water molecule using gas-phase geometries (A) and solvent-phase (water) geometries (B).

migration of the uracil anion away from the sugar moiety with no association of the water nucleophile. It is interesting to note that the dissociative pathway on the solvent-phase optimized PES is very flat (Figure 4C), where the intermediate occurs in a mere 2 kJ mol-1 well on the surface that is too shallow to see at the 5 kJ mol-1 resolution of the figure. However, the corresponding Gibbs energy surface (Figure 6B) indicates that there is a dissociative intermediate in a well with a depth of over 15 kJ mol-1 (light blue region). This suggests that in addition to allowing for comparison to experiment, the Gibbs energy surface can be used to aid the identification of important points that may be difficult to locate on the potential energy surface. A similar conclusion was drawn in a previous study of the glycosidic bond dissociation in 2′deoxyguanosine due to platination.24 It is also of note that the intermediate identified on the potential energy surface appears at a nucleophile distance of 2.6 Å, while the intermediate on the Gibbs energy surface occurs at a larger nucleophile distance (3.2 Å), which implies less association of the nucleophile. Thus, consideration of the Gibbs surface can also change the structure of intermediates along the reaction coordinate. Contrary to the gas-phase model where the SN2 (red) region extends to a nucleophile distance of ∼2.8 Å (Figure 4A), the SN2 region only reaches a nucleophile distance of ∼2.4 Å in bulk water and therefore a majority of the solvent surface (Figure

Przybylski and Wetmore 4C) is indicative of an SN1 mechanism. Nevertheless, a concerted, SN2, transition state can be identified in the solvent phase (TSC, Figure 7B), which leads to an approximate activation energy of 175 kJ mol-1. Compared with the gas-phase structure (TSC, Figure 5B), this SN2 transition state is less strained, where the strong Ow-H · · · O2 hydrogen bond seen in the gas phase is not present in the solvent-phase optimized geometry. Instead, the water nucleophile initiates attack on C1′ from directly below the sugar moiety with no explicit activation by uracil, and the uracil anion begins to depart from the cationic sugar with no explicit charge stabilization. Similar to the gas-phase product (P, Figure 5B), the solventphase product (P, Figure 7) shows transfer of a proton from the nucleophilic water to O2 of the uracil anion, where the O2 protonated uracil tautomer then interacts with the C1′-hydroxyl group through an O2-H · · · Ow hydrogen bond. However, unlike the gas-phase optimized product, this structure is not connected to either the SN1 or the SN2 transition states on the solventphase surface. In the case of the SN1 pathway, structures with nucleophile distances