Ab Initio ONIOM−Molecular Dynamics (MD) Study on the Deamination

We applied the ONIOM−molecular dynamics (MD) method to the hydrolytic deamination of cytidine by cytidine deaminase, which is an essential step of t...
0 downloads 0 Views 730KB Size
Download PDF
J. Phys. Chem. B 2007, 111, 9965-9974

9965

Ab Initio ONIOM-Molecular Dynamics (MD) Study on the Deamination Reaction by Cytidine Deaminase Toshiaki Matsubara,*,† Michel Dupuis,‡ and Misako Aida† Center for Quantum Life Sciences and Graduate School of Science, Hiroshima UniVersity, 1-3-1, Kagamiyama, Higashi-Hiroshima 739-8530, Japan, Pacific Northwest National Laboratory, Chemical Sciences DiVision, K1-83, Battelle BouleVard, Richland, Washington 99352 ReceiVed: April 7, 2007; In Final Form: June 5, 2007

We applied the ONIOM-molecular dynamics (MD) method to the hydrolytic deamination of cytidine by cytidine deaminase, which is an essential step of the activation process of the anticancer drug inside the human body. The direct MD simulations were performed for the realistic model of cytidine deaminase by calculating the energy and its gradient by the ab initio ONIOM method on the fly. The ONIOM-MD calculations including the thermal motion show that the neighboring amino acid residue is an important factor of the environmental effects and significantly affects not only the geometry and energy of the substrate trapped in the pocket of the active site but also the elementary step of the catalytic reaction. We successfully simulate the second half of the catalytic cycle, which has been considered to involve the rate-determining step, and reveal that the rate-determining step is the release of the NH3 molecule.

1. Introduction The quantum mechanics/molecular mechanics (QM/MM) hybrid methodology is now routinely used for the examination of the large molecular system. When our target is a chemical reaction, the QM method that can treat the electronic state of the molecule is indispensable. However, due to the limitation of the computational time, we are very often forced to reduce the real molecule to the model molecule, which is restricted to the active region, by disregarding the environmental effects. To resolve this problem, the QM+MM hybrid method has been developed because the MM method is good enough to describe the environmental effects. The ONIOM method1-6 based on this concept, which is developed by the Morokuma group, can easily include the environmental effects into the active (QM) part by a extrapolation procedure due to the simple definition of the QM and MM layers. In contrast, the traditional QM/MM method that adopts the sophisticated link between the QM and MM regions makes the handling difficult. In this regard, the ONIOM method is more flexible and versatile compared to the conventional QM/ MM method. Therefore, the ONIOM method is now adopted in many areas of chemistry.7 Using the ONIOM method, we have revealed the environmental effects that control the reaction in the various cases of the organometallic compounds.3,8-10 In these cases, the steric effect of the substituents on the active site plays an crucial role to determine the reactivity and selectivity of the reaction. We can also utilize the ONIOM method to examine the environmental effects in the biomolecular system in the same manner. However, when our target is a biomolecule such as an enzyme, in which the property of the entire system would significantly reflect its dynamical behavior, we cannot evaluate the environmental effects exactly without including the thermal motion. * To whom correspondence should be addressed. E-mail address: [email protected]. † Hiroshima University. ‡ Pacific Northwest National Laboratory.

To examine the environmental effects of the biomolecular system by appropriately taking account of the thermal motion of the environment, we recently developed the ONIOMmolecular dynamics (MD) method by coupling the ONIOM method with the MD method.11 We have already applied the ONIOM-MD method to the cytosine and cytidine deaminase that combine with the substrates in the previous studies,11,12 and revealed that the thermal motion of the environment is an important factor of the environmental effects. The ONIOMMD simulations showed that the energy and geometry of the substrate trapped inside the pocket of the active site is strongly perturbed by the steric effect of the amino acid residues, which sandwich the substrate, due to their thermal motion. Cytidine deaminase as well as cytosine deaminase has been examined in detail from a pharmaceutical point of view because it plays the most important role in the activation process of the anticancer drug. Recently, oral therapy has been popularized for metastatic breast and colorectal cancer. The administered prodrug undergoes some conversion processes by enzymes step by step inside the human body and is finally converted to 5-fluorouracil (5-FU) that attacks tumor tissues when it reaches the cancer cell. Cytidine deaminase catalyzes the deamination of the intermediate 5′-deoxy-5-fluorocytidine (5′-DFCR) and produces 5′-deoxy-5-fluorouridine (5′-DFUR) that is the precursor of 5-FU. This is the essential step of the entire process. Cytidine deaminase is also an excellent catalyst for the deamination of cytidine.13,14 The reaction mechanism in the active site proposed on the basis of the experimental results is presented in Figure 1. Uridine and ammonia is produced from cytidine and water, where it is considered that the reaction undergoes an intermediate B, usually referred to as the tetrahedral intermediate.14 The participation of Glu104 in the deamination reaction has been speculated from the observed crystal structures that combine with the substrates.15,16 The incoming H2O molecule is first trapped at the empty site of the Zn atom in E. One of the H is abstracted from the coordinated H2O by the glutamic acid and is transferred to the cytidine,

10.1021/jp072732k CCC: $37.00 © 2007 American Chemical Society Published on Web 07/28/2007

9966 J. Phys. Chem. B, Vol. 111, No. 33, 2007

Matsubara et al. method is discussed in Section 3.3. In Section 3.4, the final step, the release of NH3 is discussed on the basis of the results calculated by the QM-MD method for the small model of the active site used in Section 3.1 to give a further insight into the rate-determining step. Conclusions are summarized in the final section. 2. Computational Details

Figure 1. Proposed catalytic cycle of the hydrolytic deamination of cytidine in the active site of cytidine deaminase.

which induces the nucleophilic attack of the hydroxyl anion on the Zn atom to the cytidine (E f A). As a result, the tetrahedral intermediate B is formed. This is the first process of the catalytic cycle. In the second process of the catalytic cycle, the -OH hydrogen in B is carried toward the -NH2 group by the glutamic acid to form NH3 and the product uridine (B f C f D). This mediation of the glutamic acid is thought to lower the energy barrier of deamination reaction. We have also recently shown that the glutamic acid significantly contributes to lower the energy barrier of the deamination reaction by means of the density functional theory (DFT) method using a small model of the active site,17 where the environmental effects of the amino acid residues in the pocket of the active site are ignored. The previous experimental studies have indicated that the second process of the catalytic cycle after the formation of the tetrahedral intermediate B is the rate-determining step for both enzymatic deamination of cytidine and cytosine.14,18 Although the reaction mechanism is being clarified, there still remain a lot of questions to understand why cytidine deaminase shows such high catalytic activity for the deamination. From the viewpoint of improving the efficacy of the anticancer drug, further insights into the function and catalysis of the enzyme are strongly desired. The contribution of the thermal motion of the neighboring amino acid residues environment to the catalysis is now one of subjects of interest. The MD study on the catalytic reaction would bring us a new finding and an understanding concerning the enzymatic reaction. In the present study, we followed the second half of the catalytic reaction, which is considered to involve the rate-determining step, for a realistic model of the enzyme using the ONIOM-MD method and including the environment and its thermal motion. Following the explanation of the computational details in Section 2, at first the potential energy surface and the intermediates and the transition states of the second half of the hydrolytic deamination of cytidine by cytidine deaminase calculated by the quantum mechanical (QM) method using the small model of the active site are discussed in Section 3.1. In the subsequent Section 3.2, we will discuss the optimized structure of the realistic model of the tetrahedral intermediate 12, which is the starting enzyme-substrate complex of the second process of the deamination. The simulation of the deamination reaction of the realistic model of cytidine deaminase by the ONIOM-MD

2.1. Realistic Model of the Enzyme and the Small Model of the Active Site. We used the crystal structure of cytidine deaminase that binds with the substrate uridine, which are registered in the protein data bank (PDB) with the ID code of 1AF2,16 as the realistic model of the enzyme for the ONIOMmolecular dynamics (MD) simulations. The tetrahedral intermediate 12 of the realistic model, from which the ONIOMMD simulations were started, was constructed based on the crystal structure of 1AF2 modifying the reaction region. For the QM part of the two-layered ONIOM, we cut out the active region from the active site of the real enzyme and formed the model system, CH3COO- + H2O + (cytosine)Zn(SH)2(H2O), capping the dangling bonds with the H atoms. We included in the model system the -COO- group of Glu104 that mediates the deamination reaction and the three-coordinate Zn complex that binds with the substrate besides the substrate and the H2O molecule. For the substrate in the model system, the ribose ring of cytidine was replaced by the H atom and cytosine was used instead of cytidine. To mimic the lone electron pair of the N of the histidine that coordinates to the Zn atom in the real system, a H2O molecule was used in the model system. For the quantum mechanical (QM) calculations to determine the potential energy surface of the deamination reaction and the structure of the intermediates and the transition states, we adopted the small model, cytosine + CH3COO- + H2O + Zn, which was used in the previous study by the density functional theory (DFT) method.17 The same small model of the active site, cytosine + CH3COO- + H2O + Zn, was also used for the QM-MD simulations. 2.2. Quantum Mechanical (QM) Calculation. The quantum mechanical (QM) calculations were performed using the GAUSSIAN03 program19 for the small model of the active site, cytosine + CH3COO- + H2O + Zn, at the Hartree-Fock (HF) level of theory with the 3-21G basis set for all the atoms to examine the shape of the potential energy surface of the deamination reaction and the features of the geometry of the active part of the intermediates and the transition states involved in the reaction. Here, we confirmed that the HF level is good enough to describe the potential energy surface and the structure of the intermediates and the transition states qualitatively by comparing to those calculated at the B3LYP level of theory in the previous study.17 The atomic charge was obtained by the natural bond orbital (NBO) analysis.20 2.3. ONIOM-Molecular Dynamics (MD) Simulation. We used the ONIOM-molecular dynamics (MD) method implemented in the HONDO-2001 program.11 In the ONIOM-MD method, the direct MD simulations were performed with the Beeman algorithm21 by calculating the ONIOM energy and its gradient on the fly. The simulations were run under the constant temperature through the use of Berendsen’s velocity scaling algorithm22 with a time step of 1 fs. We used the two-layered ONIOM methodology, where the entire system was divided into two layers of the core and outer parts as mentioned in Section 2.1. The QM calculation of the core part (the model system) were performed at the HF level of theory with the basis set, BSI, which consists of a double-ζ (14s,8p,5d)/[5s,2p,2d] set for

ONIOM-MD on Deamination by Cytidine Deaminase

J. Phys. Chem. B, Vol. 111, No. 33, 2007 9967

Figure 2. Change in the temperature (A) and the ONIOM energy (B) of the realistic model of cytidine deaminase combined with the substrate in the water solvent in the ONIOM-MD simulation at 298.15 K.

the Zn atom23 and 3-21G for the other atoms. For the MM calculation of the outer part, i.e., the other part of the enzyme and the water solvent, the AMBER9924-28 and TIP3P29 force field parameters were used. About 2500 water molecules were randomly placed around the enzyme to mimic the water solvent, and no periodic or constrained spherical boundary conditions were used. The motion of the C6-O5-Zn(S-)2(H2O) of the QM part was fixed in the MD simulations. The conformation of the enzyme was monitored during the MD simulations. The potential energy of the entire system (ONIOM energy) as well as the temperature became almost constant after 5 ps, as presented in Figure 2, and no significant conformation change of the enzyme was found during the simulation, which is shown explicitly by the change in the potential energy of the entire system. The shape of the pocket of the active site was also not changed, as shown by its mean square displacement (MSD) and the H-bonds between the substrate and the neighboring amino acid residues that were maintained during the simulation (see the Results and Discussion). The average and standard deviation of the geometric parameters and the QM energy of the substrate in the ONIOMMD simulations were calculated from the data collected every 10 fs. 2.4. QM-Molecular Dynamics (MD) Simulation. The direct MD simulations were also performed for the small model of the active site, cytosine + CH3COO- + H2O + Zn, by calculating the energy and its gradient by the QM method on the fly. The time evolution of the nuclei was performed using the Beeman algorithm21 with a time step of 1 fs under the constant temperature through the use of Berendsen’s velocity scaling algorithm.22 The QM calculations were performed at the HF level of theory with the basis set 3-21G for all the atoms. The motion of the Zn and the CH3 carbon of CH3COO- were fixed during the MD simulations to avoid the parting of the CH3COO- fragment from the substrate. 3. Results and Discussion 3.1. Potential Energy Surface of the Deamination Reaction for the Small Model of the Active Site. We at first examined the potential energy surface of the deamination reaction by the quantum mechanical (QM) method using the small model of the active site, cytosine + CH3COO- + H2O + Zn, which was adopted in the previous study.17 The potential energy surface

of the deamination reaction calculated at the HF/3-21G level is presented in Figure 3, together with the optimized structures of the intermediates and the transition states involved in the reaction. The shape of the potential energy surface of the reaction as well as the structure of the intermediates and the transition states qualitatively reproduces that calculated previously at the B3LYP/6-31G** level of theory, which suggests that the HF/3-21G level is good enough to calculate the structure and the energy of the active part. In the second half the catalytic cycle, which is considered to be the rate-determining step, uridine and ammonia are produced by the migration of the -OH hydrogen to the -NH2 group as a proton. Although this process itself has a high-energy barrier, its energy barrier is significantly lowered by the mediation of the glutamic acid.17 The reaction starts from the tetrahedral intermediate 12, where the Zn atom binds with the O5 and the -COO- group of the glutamic acid bridges the N7H13 and the O5H25 of the substrate. The H25 is abstracted as a proton by the -COO- group of the glutamic acid in the step, 12 f 13. The abstracted H25 in 13 is carried to the -NH2 group through the transition state TS5 to from the intermediate 14. In 14, the H25 has an electrostatic interaction with the N22 atom with the short distance of 1.597 Å. The H25 is transferred to the -NH2 group as a proton in the next step, 14 f 15, through the transition state TS6, and uracil and NH3 are formed in 15. Thus, the glutamic acid mediates the proton migration. The H-bond O16‚‚‚H13-N7 was maintained throughout the reaction, although the H13 migrated to the O16 in 15. The highest point of the potential energy surface of the deamination reaction was at TS6, which indicates that the formation and release of the NH3 is the rate-determining step. We have found that the rotation of the -NH2 group is a driving force of the proton migration in the step from 12 to 13.17 We examined this step in more detail in the present study. The change in the O5-H25 and H25‚‚‚O18 distances by the rotation of the -NH2 group is presented in Figure 4B. Here, the -NH2 group was rotated counterclockwise starting from 12. When the -NH2 group was rotated by 75°, the H25 atom migrated from the O5 atom to the O18 atom. When the -NH2 group was further rotated by 225°, the H25 atom came back to the O5 atom from the O18 atom. One will notice that the H25 atom migrates to the O18 atom when the lone electron pair of the -NH2 group was directed toward the O18 atom that has a large negative charge. When the lone electron pair of the -NH2

9968 J. Phys. Chem. B, Vol. 111, No. 33, 2007

Matsubara et al.

Figure 3. Potential energy profile (kcal/mol) and the optimized structures (Å) of the intermediates and the transition states for the second half of the catalytic cycle of the hydrolytic deamination of the model substrate cytosine in the model active site of cytidine deaminase at the HF/3-21G (normal type) and B3LYP/6-31G** (italic type)17 levels. The imaginary frequencies (cm-1) for the transition states are also presented together.

Figure 4. Change in the potential energy (A) and the distances of the H-bond O5-H25‚‚‚O18 and the NBO charge of the -NH2 group (B) by the counterclockwise rotation of the -NH2 group in 12. In (A), the solid and dotted lines are at the HF/3-21G and B3LYP/6-31G**17 levels, respectively.

group was directed toward the opposite direction, the H25 atom returns to the O5 atom. Therefore, it would be rational to think that the O18 abstracts the H25 as a proton to avoid the strong electron repulsion between the -COO- and the lone electron pair of the -NH2. In fact, the NBO charge of the -NH2 group was sensitive to its rotation and the negative charge increased while the lone electron pair of the -NH2 group was directed toward the migrated proton (H25), as shown in Figure 4B. The change in the potential energy by the counterclockwise rotation of the -NH2 group is presented in Figure 4A. When the -NH2 group was rotated by around 150°, the intermediate 13 is formed. This potential energy surface shows that 13 is formed from 12 by not the counterclockwise rotation but the clockwise rotation of the -NH2 group that undergoes the energetically lower transition state. 3.2. Optimized Structure of the Tetrahedral Intermediate 12 for the Realistic Model. We fully optimized the tetrahedral intermediate 12 for the realistic model of cytidine deaminase without any constraint except for the Zn(S)2(NCC) part. The optimization was performed in the water solvent at the ONIOM-

(HF/BSI:AMBER99+TIP3P) level of theory by solving the Newton’s equations of motion at 0 K. The initial geometry was prepared from the crystal structure of cytidine deaminase combined with the substrate uridine, which is registered in the PDB with the ID code of 1AF2. Only the substrate trapped in the active site, and the neighboring amino acid residues of the optimized structure of the realistic model are displayed in Figure 5. The tetrahedral intermediate was successfully reproduced for the realistic model. The Zn atom binds with the O5 atom with the distance of 1.985 Å, and the -COO- of Glu104 interacts with the H25 and H13 atoms with the distances of 1.506 and 1.703 Å, respectively. The -NH2 group is rotated clockwise by 140° compared to that in the case of the small model, which suggests that the rotation of the -NH2 group is strongly affected by the environmental effects. It should be noted here that the main chain of the CO oxygen of Tyr126 interacts with the H23 with the distance of 2.829 Å. Three H-bonds between the substrate and the amino acid residues, Ala103, Asn89, and Glu91, found in the crystal structure, were also successfully reproduced. The carbonyl group of the six-membered ring forms

ONIOM-MD on Deamination by Cytidine Deaminase

J. Phys. Chem. B, Vol. 111, No. 33, 2007 9969

Figure 5. Structure and its illustration of the tetrahedral intermediate 12 for the realistic model of cytidine deaminase in the water solvent optimized by the ONIOM-MD method. Only the active site that combines with the substrate is displayed. The other amino acid residues and the water solvent are omitted for clarity. The optimized geometric parameters (Å) are also presented together.

TABLE 1: Occurrence (%) of the H-bonds between the Substrate and the Neighboring Amino Acid Residues and the Average and Standard Deviation of the Geometric Parameters (in Å and deg) of Their H-bonds in the ONIOM-MD Simulation of the Realistic Model at 298.15 Ka ∠(O‚‚‚H-X)

d(O‚‚‚H) amino acid residue occurrence Ala103 Asn89 Glu91

76.3 46.9 100.0

average

standard deviation

average

standard deviation

1.925(1.842) 2.057(1.893) 1.535(1.465)

0.113 0.207 0.088

159.6(140.7) 163.9(162.2) 170.1(173.0)

8.9 8.9 5.2

a

These values are calculated using the geometric parameters of the H-bonds from 5 to 20 ps. The values in parentheses are the optimized geometric parameters. See Figure 5 for the H-bonds.

the H-bond with Ala103. One of the -OH groups attached to the ribose ring forms the H-bonds with Asn89 and Glu91. Here, Glu91 more strongly interacts with the substrate through the H-bond than Ala103 and Asn89, as shown by the short O‚‚‚H distance of 1.465 Å. These H-bonds prevent the substrate from relaxing when the substrate is sterically perturbed by other amino acid residues.11 Thus, these H-bonds play an important role to destabilize the substrate in energy. 3.3. ONIOM-Molecular Dynamics (MD) Simulation of the Realistic Model of the Enzyme. The ab initio ONIOMMD simulations of the realistic model of cytidine deaminase were performed at 298.15 K starting from the optimized structure of the tetrahedral intermediate 12 mentioned in the previous section. Three H-bonds between the substrate and the amino acid residues, Ala103, Asn89, and Glu91, shown in Figure 5, were maintained during the MD simulation, although the O‚‚‚H distances were lengthened by 0.070-0.164 Å compared to those of the optimized structure (Table 1). The tendency in the distance and angle of three H-bonds found in the optimized structure did not change in the ONIOM-MD simulation. We calculated the occurrence of three H-bonds on the basis of the geometries during the MD simulation with the requirements of the O‚‚‚H distance less than 2.0 Å and the ∠O‚‚‚H-X angle more than 120°. The occurrence of three H-bonds decreases in the order, Glu91 > Ala103 > Asn89, as expected from the average of the O‚‚‚H distances, which reflects

the strength of the H-bond. The calculated occurrence shows that Ala103 as well as Glu91 strongly interacts with the substrate through the H-bond, whereas Asn89 weakly interacts. In fact, Asn89 comes close to and goes away from the substrate repeatedly due to its thermal fluctuation. The mean square displacements were calculated for the local units of the substrate trapped in the active site and the neighboring amino acid residues as follows.

DMSD )

1

N

∑ |Ri(t) - Ri(0)|2 3N i

(1)

Here, N is the number of the atoms included in each unit. As presented in Figure 6A, the thermal motion of Tyr126 and Phe71, which are outside of the pocket of the active site and contact with the water solvent, is remarkably large. Both sixmembered rings of Tyr126 and Phe71 were not rotating freely but kept fluctuating within the range of the rotation angle from -11 to 11°. This large thermal motion of Tyr126 and Phe71, which sandwich the substrate with His102 (see Figure 5), would affect the substrate fixed by the H-bonds with other amino acid residues through the steric contact as we previously reported.11 In fact, the standard deviation of the QM energy of the substrate in the MD simulation was about two times larger for the realistic model (7.4 kcal/mol) than for the small model without the amino acid residues environment (4.2 kcal/mol). The mutual steric contact between Tyr126 and Phe71 would also make the displacement of both Tyr126 and Phe71 large. As presented in Figure 7A, the intermediate 12 was converted to the intermediate 13 at 13.5 ps. The -NH2 group rotates clockwise by taking the energetically favorable channel as mentioned in Section 3.1 (see Figure 4A). When the lone electron pair of the -NH2 is directed toward the -COO- of Glu104, the H25 migrates from the O5 to the O18. Although the magnitude of the oscillation of the O18‚‚‚H25 distance is large, the H25 never migrates to the O18 as long as the lone electron pair of the -NH2 is not directed toward the -COOof Glu104 by the rotation of the -NH2. The migration of the H25 takes place at the same time as the rotation of the -NH2, as shown by the change in the magnitude of the geometric

9970 J. Phys. Chem. B, Vol. 111, No. 33, 2007

Matsubara et al.

Figure 6. Change in the mean square displacement of the local units, the substrate and the neighboring amino acid residues, (A), the distances between the -NH2 group and Tyr126 (B), and the kinetic energy of the H atoms of the -NH2 group (C) in the ONIOM-MD simulation of the realistic model at 298.15 K. In (A), green: Asn89; blue: Phe71; red: ribose ring of the substrate; black: Glu91; pink: His102; pale blue: Ala103; yellow: Glu104; orange: Tyr126. In (B), red: d1; green: d2; blue: d3. See Figure 5 for the distances, d1, d2, and d3. In (C), blue: H23; red: H24.

Figure 7. Change in the distances of the H-bond O5-H25‚‚‚O18 and the dihedral angle ∠H24-N22-C6-O5 in the ONIOM-MD simulations of the realistic model at 298.15 K (A) and 400 K (B). The color of the lines represents the geometric parameters as follows; blue: d(O5-H25); red: d(H25‚‚‚O18); green: ∠H24-N22-C6-O5.

parameters at 13.5 ps. The kinetic energy of both H23 and H24 of the -NH2 increased during 1.8 ps just before the rotation of the -NH2 (from 11.6 to 13.4 ps), as presented in Figure 6C. In contrast to the case of the realistic model, the -NH2 group freely rotates around the C6-N22 axis in the case of the small model without the amino acid residues environment. The change in the geometric parameters in the QM-MD simulation of 12 for the small model, cytosine + CH3COO- + H2O + Zn, at 298.15 K, is presented in Figure 8, comparing to that in the ONIOM-MD simulation of 12 for the realistic model at 298.15 K. The frequent switch of the magnitude of the H24‚‚‚O5 (blue line) and H23‚‚‚O5 (red line) distances in Figure 8B shows that the -NH2 group rotates round many times during the MD simulation. Consequently, the H25 migrates between the O18 and O5 atoms again and again, as shown by the vigorous change in the O5-H25 (blue line) and O18-H25 (red line) distances in Figure 8A. The change in the same parameters, H24‚‚‚O5 and H23‚‚‚O5 distances, for the realistic model, is presented in Figure 8D for comparison. The -NH2 group stays at a certain position without the significant rotation before and after 13.5 ps in the case of the realistic model. It is obvious that the rotation of the -NH2 is restricted by the amino acid residues environment. To give a deeper insight into this point, we examined the interaction between the -NH2 and the neighboring amino acid residues. As a result, it was revealed that the rotation of the -NH2 group

is locked by the attractive interaction of the H with the CO oxygen of the main chain of Tyr126. The magnitude of the oscillation of the distances, d2 and d3, in Figure 5, increases at 14-15 ps just after the rotation of the -NH2, as shown in Figure 6B, because they are affected by the motion of the -NH2. On the other hand, one of the H of the hydroxyphenyl group of Tyr126 keeps staying in the vicinity of the -NH2 with the d1 distance of less than 2.5 Å, which would facilitate the rotation of the -NH2 through the steric contact due to the thermal motion. The O16‚‚‚H13-N7 H-bond between the -COO- of Glu104 and the substrate was maintained in the realistic model, as shown in Figure 8E, but not in the small model, as shown in Figure 8C. The formation and dissociation of the O16‚‚‚H13-N7 and O16‚‚‚H23-N22 H-bonds were repeated alternately in the small model due to the large thermal motion of the O16 atom. On the other hand, the position of the O16 and the O18 was kept so as to always form the O16‚‚‚H13-N7 and O18‚‚‚H25-O5 H-bonds in the realistic model. As a matter of course, the motion of the -COO- of Glu104 is restricted because the space inside the pocket of the active site is limited. However, the space provided for Glu104 seems to be well designed to permit the O18 to merely shuttle between the O5 and the N22 as a carrier of the proton. This would be responsible for the fact that Glu104 efficiently functions as a mediator of the proton transfer.

ONIOM-MD on Deamination by Cytidine Deaminase

J. Phys. Chem. B, Vol. 111, No. 33, 2007 9971

Figure 8. Change in the selected distances in the QM-MD simulation of the small model (A,B,C) and in the ONIOM-MD simulation of the realistic model (D,E) at 298.15 K. The color of the lines represents the distances as follows; in (A), blue: d(O5-H25); red: d(O18-H25), in (B,D), blue: d(O5-H24); red: d(O5-H23), and in (C,E), blue: d(N7-O16); red: d(N22-O16).

At 298.15 K, only the migration of the H25 from the O5 to the O18 (12 f 13) proceeded by the rotation of the -NH2 group, and the subsequent reaction never proceeded during the MD simulations. This would be attributed to the fact that the transition states in the subsequent steps, TS5 and TS6, lie higher in energy than the transition state of the step, 12 f 13. The thermal motion of the active part would not be enough to readily go over the energy barrier of the subsequent steps. According to the potential energy surface of the deamination reaction calculated for the small model by the quantum mechanical method (Figure 3), at least, the energy of 10.2 kcal/mol is required to pass through the transition state TS6, which is at the highest point. On the other hand, in the MD simulation, the magnitude of the oscillation of the potential energy of the QM part of the realistic model is 14.8 kcal/mol. This suggests that the deamination reaction will be in principle completed when more than 70% of the energy of the QM part concentrates on the normal mode corresponding to the reaction coordinate. However, to complete the deamination reaction within the limited time of the MD simulation, we tried to increase the probability of the occurrence of the deamination reaction raising the temperature of the system. When the temperature of the entire system is raised up to 400 K, the intermediate 12 is transformed to the intermediate 13 by the rotation of the -NH2 group in the earlier stage (at 4.5 ps) as presented in Figure 7B. Because the rotation of the -NH2 group is locked by the attractive interaction of the H with the CO oxygen of the main chain of Tyr126 as mentioned above, the intermediate 13 does not return to the intermediate

12 easily even at 400 K. The increase in the concentration of 13 would promote the catalytic reaction. As far as the rotation of the -NH2 group is not unlocked, the transformation between 12 and 13 does not occur. The magnitude of the oscillation of the potential energy of the QM part became 1.5 times larger by the increase in the temperature of the entire system. Nevertheless, the subsequent reaction 14 f 15 did not occur. The existence ratio of the intermediate 14 was only 0.1%. Therefore, we further raised the temperature of the entire system up to 600 K to increase the probability of the formation of the intermediate 14. The transformation between 12 and 13 occurs frequently at 600 K, as presented in Figure 9A, because the -NH2 group rotates freely. The intermediate 14 was formed for a while after 9.5 ps, which was included as much as 31.2% during 3 ps from 9.6 to 12.6 ps (Table 2). We found that, besides 14, another intermediate 14′ that exists between 14 and the transition state TS6 is included at 0.3%. In the intermediate 14′, the C6-N22 bond is hardly lengthened, although the H25 is almost transferred to the N22. We have already shown in the previous study by the MO method that the intermediate 14′ can exist under certain circumstances.17 However, the intermediate 14′ is not found by the MO method for the model system cytosine + CH3COO- + H2O + Zn. After 12.6 ps, the intermediate 14 went back to the intermediate 13 without going to 15. Because further increase in the temperature of the entire system induces the conformation change of the enzyme, we instead increased intentionally the kinetic energy of the local unit that includes the reaction part using the Berendsen’s velocity scaling

9972 J. Phys. Chem. B, Vol. 111, No. 33, 2007

Matsubara et al.

Figure 9. Change in the selected distances of the QM part in the ONIOM-MD simulations of the realistic model at 600 K. In (B), the kinetic energy of the local unit consisting of the N7, C8, C9, N10, C11, O12, N22, H23, and H24 atoms was intentionally increased after 9.6 ps. See the text for the details. The color of the lines represents the distances as follows; in (A), green: d(N22-H25); blue: d(O5-H25); red: d(O18-H25); black: d(O18-N22), and in (B), green: d(N22-H25); blue: d(O5-H25); red: d(O18-H25); black: d(C6-N22).

TABLE 2: Existence Ratio (%) of the Intermediates, 14 and 14′, in the ONIOM-MD Simulations of the Realistic Model at 600 Ka intermediate

case A

case B

14 14′

31.2 0.3

22.3 8.7

a These values are calculated from the data during 3 ps from 9.6 to 12.6 ps. In case B, the kinetic energy of the local unit consisting of the N7, C8, C9, N10, C11, O12, N22, H23, and H24 atoms was intentionally increased after 9.6 ps. See the text for the details.

algorithm.22 We selected nine atoms of the substrate composing the local unit, N7, C8, C9, N10, C11, O12, N22, H23, and H24. This selection is reasonable, because the thermal energy would be localized on a vibrational mode of the substrate when the reaction proceeds. The H25 and H13 atoms are already “thermally activated” at 600 K. We increased the kinetic energy of the local unit consisting of these nine atoms by two times, when 9.6 ps passed and the intermediate 14 was formed. As a result, the formed NH3 molecule on the C6 is released 3 ps later and the deamination reaction was successfully completed, as shown in Figure 9B. This result suggests that the release of the NH3 molecule from the C6 (14 f 15) is the rate-determining step. During 3 ps, from 9.6 to 12.6 ps before releasing the NH3, the intermediate 14 was included at 22.3%, and the intermediate 14′ was also included at 8.7%. Although the intermediate 14 once retuned to the intermediate 13 after increasing the kinetic energy of the local unit, 14 was formed again at 11.3 ps. The H25 showed a reciprocating motion between the O18 and the N22 during 1.4 ps (from 11.3 to 12.7 ps) just before the release of the NH3 molecule. This suggests that the NH3 molecule is released when the motions of both H25 and N22 match to each other. 3.4. QM-Molecular Dynamics (MD) Simulation of the Release of the NH3 Molecule. We examined in detail the ratedetermining step, the release of the NH3 molecule (14 f 15), by the QM-molecular dynamics (MD) method using the small model, cytosine + CH3COO- + H2O + Zn. As mentioned in the previous section for the realistic model, the MD simulation at 298.15 K starting from the intermediate 12 never produces the final products, NH3 and uracil, similarly also for the small model. Only the transformation between the intermediates 12 and 13 occurs. When the temperature of the entire system is increased, the weak interaction between the -COO- and the

substrate is broken and the -COO- goes away from the substrate. The MD simulations starting from the intermediate 14 also gave the same results. The intermediate 14 goes to 13 or 12 without producing 15. However, we already know that the reaction would proceed when a thermal energy is provided for an appropriate local unit as shown for the realistic model in the previous section. We therefore selected three atoms, H25, H13, and N22, which would significantly contribute to the normal mode corresponding to the reaction coordinate from the formed 12 (or 13) to 15 and increased the kinetic energy of these three atoms by 2.4 times in the same manner adopted for the realistic model in the previous section when the first 5 ps passed starting from the intermediate 14 at 298.15 K. The thermal energy would be efficiently provided for the normal mode that corresponds to the reaction coordinate from 12 (or 13) to 15 by this increase in the kinetic energy of the selected three atoms, and then the probability of the occurrence of the reaction would be enhanced. As presented in Figure 10B, the intermediate 14 once returned to the intermediate 13 after 9 ps. However, the intermediate 13 is transformed to the intermediate 14 again after 14 ps, and the NH3 molecule is finally released after 18.5 ps, as expected. The existence ratio of the intermediate 14 during 12 ps from 5 to 17 ps was 46.4%, as shown in Table 3. It was also found that the intermediate 14′ exists by 1.2% without the amino acid residues environment. On the other hand, without the increase in kinetic energy of the selected three atoms, the intermediate 14 returned to the intermediate 12 or 13, as presented in Figure 10A, although the intermediate 14 was trapped during the first 6 ps. The NH3 molecule was never released. The existence ratio of the intermediate 14 during 12 ps from 5 to 17 ps decreased to 6.9% and the intermediate 14′ was not included. Thus, the QM-MD simulation of the small model also showed that the release of the NH3 molecule (14 f 15) is the rate-determining step. As far as we tried, the formed intermediate 12 or 13 from 14 did not complete the reaction, even if the kinetic energy of the -NH2 group or only the H13 and H25 atoms were intentionally increased. These results show that a localization of the thermal energy in an appropriate part would be important as a driving force of the reaction. Comparing Figure 10D with Figure 10C, it is obvious that the H-bond O16‚‚‚H13-N7 is an important factor to complete the final step, 14 f 15. The H-bond O16‚‚‚H13-N7 is maintained until the NH3 molecule is released, as shown in

ONIOM-MD on Deamination by Cytidine Deaminase

J. Phys. Chem. B, Vol. 111, No. 33, 2007 9973

Figure 10. Change in the selected distances in the QM-MD simulations of the small model at 298.15 K. In (B) and (D), the kinetic energy of the local unit consisting of the H25, H13, and N22 atoms was intentionally increased after 5 ps. See the text for the details. The color of the lines represents the distances as follows; green: d(N22-H25); blue: d(O5-H25); red: d(C6-N22); black: d(O18-H25); pink: d(H13-O16); yellow: d(N7-H13).

TABLE 3: Existence Ratio (%) of the Intermediates, 12, 13, 14, and 14′, in the QM-MD Simulations of the Small Model at 298.15 Ka intermediate

case A

case B

12 13 14 14′

5.4 39.0 6.9 0.0

3.3 21.4 46.4 1.2

a These values are calculated from the data during 12 ps from 5 to 17 ps. In case B, the kinetic energy of the local unit consisting of the H25, H13, and N22 atoms was intentionally increased after 5 ps. See the text for the details.

Figure 10D, although the magnitude of the oscillation of the O16‚‚‚H13 distance is large. On the other hand, in Figure 10C, the O16‚‚‚H13 interaction is broken, as shown by the O16‚‚‚H13 distance of more than 3 Å, and the NH3 molecule is never produced. One will notice at a glance that the magnitude of the vibration of the N7-H13 bond is explicitly large in Figure 10D compared to that in Figure 10C because the kinetic energy of the H13 is intentionally increased. The intensive stretch of the N7-H13 bond enlarges its polarization so that the electrostatic interaction of the H13 with the O16 is strengthened. This is the reason why the O16‚‚‚H13-N7 H-bond is maintained. The magnitude of the oscillation of the O16‚‚‚H13 distance is obviously smaller in 14 (from 5 to 9 ps and from 14 to 18 ps) than in 13 (from 9 to 14 ps). The H13 was abstracted as a proton by the O16 of the -COO- at the same time as the release of

the NH3 molecule, as shown by the change in the N7-H13 (yellow line) and O16‚‚‚H13 (pink line) distances, which is consistent with the result by the quantum mechanical calculation for the small model (see section 3.1). 4. Concluding Remarks We applied the ONIOM-molecular dynamics (MD) method to the hydrolytic deamination of cytidine to give a further insight into the function of cytidine deaminase by taking account of the thermal motion of the entire system. We focused on the second half of the catalytic cycle, which has been thought to involve the rate-determining step on the basis of the experimental results. The ONIOM-MD simulations were performed for the realistic model of cytidine deaminase along the reaction coordinate calculated by the quantum mechanical (QM) method for the small model of the active site. The substrate trapped in the active site is fixed by the amino acid residues inside the pocket, Ala103, Asn89, and Glu91, through the H-bonds. Other amino acid residues outside the pocket, Tyr126 and Phe71, which sandwich the substrate with His102, strongly perturbs the geometry and energy of the substrate through the steric contact due to their large thermal motion. The -COO- of Glu104 is kept at a certain position to have an electrostatic interaction with the substrate in the best way by the environment inside the pocket so that Glu104 can efficiently mediate the transfer of the proton. The abstraction of the proton by the -COO- of Glu104 in the first step, 12 f 13, is driven by the

9974 J. Phys. Chem. B, Vol. 111, No. 33, 2007 clockwise rotation of the -NH2 group in 12, and the formed intermediate 13 is stabilized in energy by the interaction with Tyr126. The ONIOM-MD simulations showed that the final step (14 f 15), the release of the NH3 molecule, is the ratedetermining step. We succeeded in simulating the rate-determining step by increasing the kinetic energy of a local unit of the reaction part of the substrate. In the rate-determining step, we found a short-lived intermediate 14′ besides 14. The present study of the realistic model of cytidine daminase by the ONIOM-MD method including the thermal motion successfully showed that the neighboring amino acid residue affects not only the geometry and energy of the substrate but also the elementary step of the catalytic cycle and is an important factor of the environmental effects to support the deamination reaction. Acknowledgment. T.M. and M.A. were supported in part by grants from the Ministry of Education, Culture, Sports, Science and Technology of Japan. M.D. was supported by the Division of Chemical Sciences, Office of Basic Energy Sciences, and by the Office of Biological and Environmental Research of the U.S. Department of Energy (DOE). Battelle operates Pacific Northwest National Laboratory for DOE. Supporting Information Available: Listings giving the optimized Cartesian coordinates of the equilibrium structures, 12, 13, 14, and 15, and the transition states, TS5 and TS6, in Figure 3 at the HF/3-21G level. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Maseras, F.; Morokuma, K. J. Comput. Chem. 1995, 16, 11701179. (2) Matsubara, T.; Sieber, S.; Morokuma, K. Int. J. Quantum Chem. 1996, 60, 1101-1109. (3) Matsubara, T.; Maseras, F.; Koga, N.; Morokuma, K. J. Phys. Chem. 1996, 100, 2573-2580. (4) Svensson, M.; Humbel, S.; Froese, R. D. J.; Matsubara, T.; Sieber, S.; Morokuma, K. J. Phys. Chem. 1996, 100, 19357-19363. (5) Dapprich, S.; Koma´romi, I.; Byun, K. S.; Morokuma, K.; Frisch, M. J. J. Mol. Struct. (THEOCHEM) 1999, 461-462, 1-21. (6) Vreven, T.; Morokuma, K. J. Comput. Chem. 2000, 21, 14191432. (7) Morokuma, K. Bull. Korean Chem. Soc. 2003, 24, 797-801. (8) Nozaki, K.; Komaki, H.; Kawashima, Y.; Hiyama, T.; Matsubara, T. J. Am. Chem. Soc. 2001, 123, 534-544. (9) Sahnoun, R.; Matsubara, T.; Yamabe, T. Organometallics 2000, 19, 5661-5670.

Matsubara et al. (10) Nozaki, K.; Sato, N.; Tonomura, Y.; Yasutomi, M.; Takaya, H.; Hiyama, T.; Matsubara, T.; Koga, N. J. Am. Chem. Soc. 1997, 119, 1277912795. (11) Matsubara, T.; Dupuis, M.; Aida, M. Chem. Phys. Lett. 2007, 437, 138-142. (12) For cytosine deaminase, see: Matsubara, T.; Dupuis, M.; Aida, M. J. Comput. Chem. 2007, in press. (13) Snider, M. J.; Gaunitz, S.; Ridgway, C.; Short, S. A.; Wolfenden, R. Biochemistry 2000, 39, 9746-9753. (14) Snider, M. J.; Reinhardt, L.; Wolfenden, R.; Cleland, W. W. Biochemistry 2002, 41, 415-421. (15) Carlow, D. C.; Smith, A. A.; Yang, C. C.; Short, S. A.; Wolfenden, R. Biochemistry 1995, 34, 4220-4224. (16) Xiang, S.; Short, S. A.; Wolfenden, R.; Carter, C. W., Jr. Biochemistry 1997, 36, 4768-4774. (17) Matsubara, T.; Ishikura, M.; Aida, M. J. Chem. Inf. Model. 2006, 46, 1276-1285. (18) Yao, L.; Li, Y.; Wu, Y.; Liu, A.; Yan, H. Biochemistry 2005, 44, 5940-5947. (19) 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.; Bakken, V.; 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.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03; Gaussian, Inc.: Pittsburgh, PA 2003. (20) Glendening, E. D.; Reed, A. E.; Carpenter, J. E.; Weinhold, F. NBO, version 3.1. (21) Beeman, D. J. Comput. Phys. 1976, 20, 130-139. (22) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684-3690. (23) Scha¨fer, A.; Horn, H.; Ahlrichs, R. J. Chem. Phys. 1992, 97, 25712577. (24) Aqvist, J. J. Phys. Chem. 1990, 94, 8021-8024. (25) Ross, W. S.; Hardin, C. C. J. Am. Chem. Soc. 1994, 116, 60706080. (26) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179-5197. (27) Moyna, G.; Williams, H. J.; Nachman, R. J.; Scott, A. I. Biopolymers 1999, 49, 403-413. (28) Wang, J.; Cieplak, P.; Kollman, P. A. J. Comput. Chem. 2000, 21, 1049-1074. (29) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926-935.