1514
2007, 111, 1514-1516 Published on Web 01/31/2007
Theoretical Examination of Two Opposite Mechanisms Proposed for Hepatitis Delta Virus Ribozyme Kai Wei, Lei Liu,* Yu-Hui Cheng, Yao Fu, and Qing-Xiang Guo* Department of Chemistry, UniVersity of Science and Technology of China, Hefei 230026, China ReceiVed: January 7, 2007
Two mechanisms were previously proposed for the hepatitis delta virus (HDV) ribozyme where an activesite cytosine residue (C75) either functioned as a general base to deprotonate the 2′-OH at the rupture site or as a general acid to protonate the O5′ leaving group. Here, we reported the first theoretical examination of the two mechanisms using a combination of the quantum mechanics (QM)/molecular mechanics (MM), molecular dynamics (MD), and near-attack-conformation (NAC) techniques. Our theoretical results supported the C75-acid mechanism, which was demonstrated to have an unfavorable starting geometry (in agreement with the crystallographic data) but a significantly lower energy barrier as compared to the C75-base mechanism. Therefore, the chemical details of the transition state in the HDV ribozyme may dramatically differ from those inferred from the structural studies.
Introduction Ribozymes are biological catalysts comprised of folded RNAs. They accelerate chemical reactions using many of the same strategies (such as general acid/base catalysis and substrate desolvation) as protein enzymes.1 Nonetheless, because RNAs do not have those well-characterized catalytic groups available to proteins, many details about how the ribozymes achieve catalysis remain enigmatic.2 For instance, in the recent studies about the hepatitis delta virus (HDV) ribozyme, an active-site cytosine residue (C75) was proposed to facilitate the proton transfer.3,4 However, it remains a matter of debate (Figure 1) whether C75 acts as a general base to deprotonate the 2′-OH at the rupture site (as suggested by the crystallographic studies)3 or as a general acid to protonate the O5′ leaving group (as suggested by reactivity-pH profiles and chemical substitution studies).4 Here, we report, for the first time, a theoretical examination of the two contradictory mechanisms proposed for HDV. The initial structure of HDV was chosen from the X-ray data reported for the C75U mutant (PDB code 1SJ3). After the U75 residue was changed back to cytosine, hydrogens of the heavy atoms were added using the Psfgen program in the NAMD software.5 The structure (with one crystallographic Mg2+ ion) was relaxed in the gas phase using the Parm99 force field.6 The relaxed structure was immersed in a cubic TIP3P water box (70 × 110 × 70 Å3) and neutralized by addition of Na+ couterions. The resulting system (total number of atoms ) 46 770) was equilibrated through a nine-step molecular dynamics protocol, and the equilibrated structure was subjected to a unconstrained molecular dynamics simulation with periodical boundary conditions (time step ) 1.0 fs, cutoff radius for nonbonding interaction ) 12 Å, isothermic-isobaric ensemble, * To whom correspondence should be addressed. E-mail:
[email protected] (L.L.);
[email protected] (Q.-X.G.).
10.1021/jp070120u CCC: $37.00
Figure 1. General acid/base catalysis by the HDV ribozyme: (a) proposed mechanism for general acid catalysis of phosphodiester cleavage by C75; (b) proposed mechanism for general base catalysis by C75.
T ) 300 K, P ) 1 atm). The electrostatic interactions were calculated by the particle mesh Ewald method.7 The SHAKE algorithm8 was used to freeze all of the bond lengths involving hydrogens. To compare the two mechanisms proposed for HDV, we defined the following chemically distinctive species as shown in Figure 2. First, protonation of either the Mg-hydroxyl or C75 gives the relaxed starting material for each mechanism (i.e., “C75-base” or “C75-acid”). The initial concentration difference between these two species can be calculated from the difference between the pKa of Mg-OH2 bound in HDV and the pKa of C75-H bound in HDV. Second, by using appropriate geometry constraints, we can define the near-attack conformation (NAC)9 for each mechanism (i.e., NAC-C75-base and NAC-C75-acid). The NAC concentrations are calculated on the basis of empirical criteria to determine the effective concentration that comes out naturally in a Potential of Mean Force simulation. The free energy of these NAC species can be calculated by counting their probability of appearance in the molecular dynamics (MD) simulation. Finally, TS-C75-base and TS-C75-acid describe the © 2007 American Chemical Society
Letters
J. Phys. Chem. B, Vol. 111, No. 7, 2007 1515
Figure 2. Chemically distinctive species in the HDV catalysis.
transition states linking the NAC and the product. The energy barrier from NAC to the TS can be estimated from quantum chemistry calculations on the active-site models. To calculate the pKa of MgOH2 bound in HDV, we used the quantum mechanics (QM)/molecular mechanics (MM) method developed by Zhang et al.10 Briefly, we defined the proton transfer reaction between two water molecules as the QM layer. Under hammonic contraints, we forced the proton to move from the Mg-coordinated water to the free water in steps, where Mg and the rest of the HDV were treated by the force field. Through the standard free energy perturbation method,10 we obtained the free energy corresponding to the proton transfer process. The pKa value for MgOH2 bound in the HDV ribozyme was thus predicted to be 9.2. Compared to the hydrolysis constant of Mg2+ in water (10-11.4), MgOH2 bound in HDV was more acidic than free MgOH22+. Using a similar procedure, we also calculated the pKa of C75-H bound in HDV to be 6.4. This value was in good agreement with the experimental value, 6.1.4 On the basis of the two pKa values, we calculated that, at pH 7.0, the concentration for the C75-base species was 0.8, whereas the concentration for C75-acid was 1.1 × 10-3 (assuming that the total concentration for HDV was 1.0). Having calculated the concentrations for C75-acid and C75base, we next utilized the NAC technique developed by Bruice et al. to estimate the effective concentrations for the near-attack conformations corresponding to the two mechanisms. In defining the NAC for the C75-base mechanism, we measured the distance (d1) between the N3 atom of the C75 residue and O2′ of the U-1 residue and the distance (d2) between the oxygen of the Mg-coordinated water and O5′ of the G1 residue. If d1 and d2 were in the range from 2.698 ( 0.4 to 2.614 ( 0.4 Å, the corresponding structure was counted as an NAC. Similarly, for the C75-acid mechanism, we measured the distance (d1′) between N3 of the C75 residue and O5′ of the G1 residue and the distance (d2′) between the oxygen atom of the Mgcoordinated OH group and O2′ of the U-1 residue. If d1′ and d2′ were in the range from 2.623 ( 0.4 to 2.543 ( 0.4 Å, the corresponding structure was counted as a NAC. By running 1 ns of unconstrainted MD for either the C75-acid or C75-base HDV ribozyme, we determined that the probability for the occurrence of NAC-C75-base was 10% relative to C75-base, whereas the probability for the occurrence of NAC-C75-acid was 0.01% relative to C75-acid. As a consequence, the effective concentration for NAC-C75-base was estimated to be 0.04, whereas the effective concentration for NAC-C75-acid was 1.1 × 10-7. Finally, on the basis of the NACs, we constructed two cutoff models for the C75-acid and C75-base mechanisms, respectively, by removing the periphery residues that were not directly interacted with C75, U-1, and the Mg-bound waters. Thus, in the reaction center, there remained the cytosine ring of C75, the U-1 residue without the nucleobase, Mg2+ and its coordina-
Figure 3. Cutoff models for the C75-base (left) and C75-acid (right) mechanisms.
Figure 4. Free energy barriers from NACs to the final products.
tion waters, the phosphate groups of C22 and U23, and the hydroxyl group of U20 (Figure 3).11 Overall, there were 42 nonhydrogen atoms in the cutoff models, each carrying a total charge of -e. The cutoff models were then subjected to full optimization using the B3LYP/3-21G(d) method. The effect of environment was modeled with the polarizable continuum solvation model by using water as the solvent.12 After geometry optimization was completed, the energy for the fully optimized structure was calculated using the B3LYP/6-31+G(d) method corrected by the zero-point energy, temperature corrections, and solvation energy. It was found that from the cutoff model to the final product the system had to undergo two transition states.13 In the first transition state, the imaginary frequency corresponds to the transfer of the proton from O2′ of the U-1 residue to the general base, coupled with the nucleophilic attack of O2′ to the scissile phosphate. This gives a trigonal bipyramidal intermediate in which the phosphorus atom is pentacoordinated. Subsequently, in the second transition state, the imaginary frequency corresponds to the P-O bond breaking event which yields the final self-cleaved product. For the C75-base mechanism, the rate determining step was found to be the second step (i.e., P-O bond cleavage) and the free energy barrier was 31.3 kcal/mol (Figure 4). For the C75-acid mechanism, the rate determining step was found to be the first step and the free energy barrier was 14.6 kcal/mol. On the basis of the concentrations for NACs and the energy barriers for the reactions from NACs to the final product, we could now calculate the overall rate constants. For the C75base mechanism (NAC concentration ) 0.08, barrier ) 31.3 kcal/mol), the rate constant was calculated to be 1.76 × 10-11 s-1 corresponding to an apparent activation energy of 32.8 kcal/ mol. For the C75-acid mechanism (NAC concentration ) 1.1 × 10-7, barrier ) 14.6 kcal/mol), the rate constant was
1516 J. Phys. Chem. B, Vol. 111, No. 7, 2007 calculated to be 2.3 × 10-5 s-1 corresponding to an apparent activation energy of 24.3 kcal/mol. Noteworthily, the experimental apparent activation energy for HDV was 19.5 kcal/mol.4 As a consequence, our theoretical results14 supported the C75acid mechanism. Noteworthily, the C75-acid mechanism has also been favored by some more recent experiments.15 Further examination of the theoretical data suggested that the NAC concentration for the C75-base mechanism (0.08) was much higher than that for the C75-acid mechanism (1.1 × 10-7). This observation was in agreement with the crystallographic study on HDV where C75 was found to stay close to 2′-OH at the rupture site in the crystal (i.e., in the starting material state).3 However, it was now clear that the C75-base mechanism must overcome a huge energy barrier of 31.3 kcal/mol, whereas the C75-acid mechanism only needs to overcome a barrier of 14.6 kcal/mol. As a consequence, the C75-acid mechanism was still the favorable one, even though it had a very unfavorable starting geometry. Thus, a valuable lesson learned from the study on HDV is that the chemical details of the transition state in the RNA catalysts may dramatically differ from those inferred from the structural studies. Acknowledgment. We thank NSFC (no. 20472079) for the financial support. Supporting Information Available: Detailed calculation procedures. This material is available free of charge via the Internet at http://pubs.acs.org.
Letters References and Notes (1) Lilley, D. M. J. Trends Biochem. Sci. 2003, 28, 495. (2) Emilsson, G. M.; Nakamura, S.; Roth, A.; Breaker, R. R. RNA 2003, 9, 907. (3) Ke, A.; Zhou, K.; Ding, F.; Cate, J. H. D.; Doudna, J. A. Nature 2004, 429, 201. (4) Nakano, S.-I.; Durge, M. C.; Philip, C. B. Science 2000, 287, 1493. (5) Canto, J.; Haro, I.; Alsina, M. A.; Perez, J. J. J. Phys. Chem. B 2003, 107, 6603. (6) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; K., S. J. Comput. Chem. 2005, 26, 1781. (7) Shields, G. C.; Laughton, C. A.; Orozco, M. J. Am. Chem. Soc. 1998, 120, 5895. (8) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089. (9) Lightstone, F. C.; Bruice, T. C. J. Am. Chem. Soc. 1994, 116, 10789. (10) Zhang, Y.; Liu, H.; Yang, W. J. Chem. Phys. 2000, 112, 3483. (11) Luptak, A. F.-D. A. R.; Zhou, K.; Zilm, K. W.; Doudna, J. A. J. Am. Chem. Soc. 2001, 123, 8447. (12) Miertus, S.; Scrocco, E.; Tomasi, J. Chem. Phys. 1981, 55, 117. (13) (a) Torres, R. A.; Himo, F.; Bruice, T. C.; Noodleman, L.; Lovell, T. J. Am. Chem. Soc. 2003, 125, 9861. (b) Gregersen, B. A.; Lopez, X.; York, D. M. J. Am. Chem. Soc. 2004, 126, 7504. (14) The difference between theoretical (24.3 kcal/mol) and experimental activation energies (19.5 kcal/mol) could originate from the following: (1) NAC simulation where the cutoff distance (i.e., 0.4 Å) was not a rigorously defined parameter or (2) energy calculation at the B3LYP/6-31G(d) level. (15) (a) Das, S. R.; Piccirilli, J. A. Nat. Chem. Biol. 2005, 1, 45. (b) Doudna, J. A.; Lorsch, J. R. Nat. Struct. Mol. Biol. 2005, 12, 395.
