Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
J Phys Chem A. Author manuscript; available in PMC 2010 November 12.
Published in final edited form as:
PMCID: PMC2783400

The hydrolysis activity of Adenosine triphosphate in myosin: a theoretical analysis of anomeric effects and the nature of transition state


Combined quantum mechanical/molecular mechanical (QM/MM) calculations with density functional theory are employed to analyze two issues related to the hydrolysis activity of Adenosine triphosphate (ATP) in myosin. First, we compare the geometrical properties and electronic structure of ATP in the open (post-rigor) and closed (pre-powerstroke) active sites of the myosin motor domain. Compared to both solution and the open active-site cases, the scissile Pγ-O bond of ATP in the closed active-site is shown to be substantially elongated. Natural Bond Orbital (NBO) analysis clearly shows that this structural feature is correlated with the stronger anomeric effects in the closed active site, which involve charge transfers from the lone-pairs in the non-bridging oxygen in the γ phosphate to the anti-bonding orbital of the scissile bond. However, an energetic analysis finds that the ATP molecule is not significantly destabilized by the Pγ-O bond elongation. Therefore, despite the notable perturbations in the geometry and electronic structure of ATP as its environment changes from solution to the hydrolysis-competent active site, ground state destabilization is unlikely to play a major role in enhancing the hydrolysis activity in myosin. Second, two-dimensional potential energy maps are used to better characterize the energetic landscape near the hydrolysis transition state. The results indicate that the transition state region is energetically flat and a range of structures representative of different mechanisms according to the classical nomenclature (e.g., “associative”, “dissociative” and “concerted”) are very close in energy. Therefore, at least in the case of ATP hydrolysis in myosin, the energetic distinction between different reaction mechanisms following the conventional nomenclature is likely small. This study highlights the importance of (i). explicitly evaluating the relevant energetic properties for determining whether a factor is essential to catalysis and (ii). broader explorations of the energy landscape beyond saddle points (even on free energy surface) for characterizing the molecular mechanism of catalysis.


The hydrolysis of nucleotides such as ATP and GTP is involved in many fundamental biological processes, especially in signal transduction and energy transduction events1. It is, therefore, useful to understand, in molecular terms, the factors that control the hydrolysis activity of nucleotides in biomolecules. Here we discuss two specific issues, which are related to the properties of ATP in the ground (prior to hydrolysis) and hydrolysis transition state, respectively. As we will show, careful theoretical calculations that explicitly probe the relevant energetics are particularly valuable in providing the key insights in both cases.

The first issue concerns whether ATP is chemically “activated” in the nucleotide binding site of an ATPase (e.g., myosin) relative to solution, making ATP more susceptible to rapid hydrolysis. This is essentially one version of the “ground state destabilization” hypothesis for enzyme catalysis, which suggests that the (electronic) structure of the substrate is perturbed significantly relative to solution by interactions with active site residues, and as a result the substrate is chemically activated for the subsequent reaction. For example, shifts of vibrational frequencies observed in infrared (IR) studies have been used to argue that the scissile P-O bond of the substrate is substantially weakened in the enzyme active site (ras-p212,3 and Ca2+-ATPase4,5) compared to the solution state; we note that for probing electronic structural changes of the substrate in an enzyme environment relative to solution, spectroscopic studies are more useful than X-ray crystallography because the latter often has to employ substrate analogues. These suggestions are in qualitative agreement with observations from several theoretical studies6-10, which reported that the optimized scissile P-O bond length (Pγ-O in ATP/GTP) in an ATPase seems to be substantially longer than the solution value (see below).

Another motivation for choosing to focus on ATP “activation” in myosin is that there is the possibility that the degree of “activation” may depend on the conformational state of the motor domain. Indeed, molecular motors11 such as myosin are striking examples where the ATP hydrolysis activity is strictly controlled to avoid futile hydrolysis, which is an important feature that governs the working efficiency of those unique systems12,13. Specifically for myosin, for example, ATP hydrolysis is believed to turn on only after the active site is switched from an open configuration to a structurally stable closed configuration14,15, and this is supported by theoretical analysis16-18. However, the precise reason for the dependence of the ATPase activity on the active site closure remains incompletely understood. It is likely that a number of factors, which include orientation of the nucleophilic water19 and electrostatic interactions20 in the active site, play a role. Whether stereoelectronic effects make a major contribution deserves a careful analysis, which requires an advanced characterization of substrate/enzyme interaction at the electronic structure level.

In this study, encouraged by our recent analysis of solute-solvent interaction for charged solute (e.g., phosphate)21, we employ the Natural Bond Orbital (NBO) approach to explore whether the ATP molecule is significantly activated in a specific conformational state of myosin. Our analyses reveal that the electronic structure and geometry of ATP undergo interesting changes among the solution state and two key conformational states of the myosin motor domain. However, we also find that the energetic variation associated with the P-O bond elongation is rather modest, thus it is unlikely that “ground state destabilization” plays a major role in the regulation of ATPase activity in myosin. This study highlights both the value and limitation of probing properties that are sensitive to electronic structure and geometry, and emphasizes the importance of directly evaluating the relevant energetic contribution for determining whether a factor is essential to catalysis.

The second issue concerns the nature of the ATP hydrolysis transition state in myosin. Conventionally, the mechanisms of phosphate hydrolysis are classified into two main extremes, the associative and dissociative mechanisms (see Scheme 1)22-24. In the associative mechanism, a nucleophilic attack by water at the phosphorous center occurs before the leaving group departs, thus a pentavalent phosphorane intermediate is involved. By contrast, in the dissociative mechanism, the leaving group dissociates from the phosphorus prior to the nucleophilic attack, consequently a metaphosphate intermediate state is involved. Under certain reaction condition, nucleophilic addition and elimination may also occur in a concerted manner, where the new P-O bond forms at the same time as the old P-O bond breaks. Specifically for ATP hydrolysis, the linear free-energy relationship (LFER) analysis of Herschlag and coworkers25 led to the conclusion that the dissociative mechanism is operational in solution, although potential ambiguity in interpreting the LFER data has also been pointed out26. For ATP hydrolysis in myosin, Trentham and coworkers studied oxygen exchange during hydrolysis using quenched-flow techniques; they concluded that the reaction goes through “a single step with direct transfer of the terminal phosphorus from ATP to water”27. More recently, by applying Raman spectroscopy to myosin·MgADP-vanadate complex, Yount and coworkers28 suggested that “this reaction proceeds close to a concerted (SN2-like) process” with an associative character. In several theoretical studies of ATP hydrolysis in myosin16,18,29, an associative mechanism is typically followed except in the recent study of Nemukhin et al.7, who advocated a rather unconventional mechanism involving Glu459 as the general base30. To help better understand the hydrolysis mechanism and the nature of transition state, we have generated a two-dimensional potential energy surface along reaction coordinates that do not specifically bias toward the associative or dissociative mechanism. We observe that the potential energy surface near the saddle point (transition state) is very flat energetically and a range of structures representative of different mechanisms according to the classical nomenclature (e.g., “associative”, “dissociative” and “concerted”) are very close in energy. Therefore, at least in the case of ATP hydrolysis in myosin, the energetic distinction between different reaction mechanisms following the conventional nomenclature is likely small.


A. ATP in the myosin active site

For the myosin system, both the post-rigor and pre-powerstroke conformational states14,15 are studied with ATP as the ligand. As discussed in previous studies14,15,31, the post-rigor state is incapable of hydrolyzing ATP, thus it was possible to solve the corresponding crystal structure with ATP as the ligand (PDB code 1FMW32,33, with the myosin motor domain from Dictyostelium discoideum). The pre-powerstroke state, on the other hand, is the conformational state in which ATP is effectively hydrolyzed, and therefore its crystal structure (PDB code 1VOM) was determined using ADPVO4 as the ligand34. Using these X-ray structures as the initial input, combined quantum mechanical/molecular mechanical (QM/MM35-37) calculations are carried out to obtain the minimized structure of the active site. Although more sophisticated sampling techniques that properly include thermal fluctuations are needed for a quantitative understanding of hydrolysis energetics (see below)18,energy minimization is sufficient for generating active site structures for the current purpose of comparing electronic structures of the substrate (ATP) in different conformational states of the motor domain.

The detailed QM/MM set-up is similar to the protocol used in our previous studies18 and therefore will not be elaborated here. The QM region (48 atoms) includes the tri-phosphate group of ATP, the Mg2+ ion and its ligands, the sidechain of the conserved Ser236 and the lytic water molecule; calculations using a larger QM region (80 atoms) that further includes residues close to the γ phosphate of ATP (Lys185, Asn233 and Ser181) have also been carried out, although the results are very similar and therefore not included. The rest of the system is treated using the CHARMM22 force field for proteins38. The structures are first minimized using the more efficient SCC-DFTBPR/MM method, where SCC-DFTBPR is an approximate density functional theory approach39 that we recently developed for treating phosphate reactions40; the low computational cost of this method helps to fully relax the system, especially the large number of MM degrees of freedom. Next, the SCC-DFTBPR/MM structures are minimized with the QM region treated at the B3LYP41-43/6-31G(d)44-47 level; the convergence threshold for the root-mean-squaregradient is set to be 0.01 kcal/(mol·Å), which is found to be sufficient for obtaining converged geometrical parameters for ATP. To analyze the electronic structure of ATP in the two myosin conformations, NBO analysis is carried out using QM/MM calculations with the QM atoms treated at the B3LYP/6-311+G* level45,48.

The SCC-DFTBPR/MM calculations are performed using a local version of CHARMM49. The B3LYP/MM calculations are carried out with the c35a1 version of CHARMM, which is interfaced with QChem (v.3.1)50,51. For NBO analysis, the NBO 5.0 package52 interfaced with QChem is used.

B. ATP in solution

Since the main goal of this work is to reveal changes in the electronic structure and energy of ATP induced by the variation of the scissile P-O bond length (Pγ-O), the solution state of ATP is studied using the tri-phsophate moiety, which is part of the QM region in the myosin set-up. To focus on effects induced by the difference in environment, the initial tri-phosphate geometry is taken from B3LYP/MM minimized structure of the myosin-ATP complex. The structure is then partially optimized in the gas-phase with the critical Pγ-O bond distance held at several values (see Sect.III for details). Single point energy and NBO analysis are then carried out at these partially optimized structures using B3LYP/6-311+G* and the PCM model for solvation53. To explore the sensitivity of the PCM results to the choice of parameters, the solvation calculations are done with both the United Atom Kohn Sham (UAKS) radii and the Pauling radii that define the solute cavity54.

C. Two-dimensional adiabatic mapping calculations

The adiabatic mapping calculations55 have only been carried out for the pre-powerstroke configuration of the active site, which has been well accepted to be the configuration where ATP is efficiently hydrolyzed15,18,31. Starting with the optimized reactant state structure obtained in Sect.II A, a two-dimensional potential energy surface is generated using adiabatic mapping. To have an unbiased description of both the associative and dissociative mechanisms, one of the reaction coordinates is defined to be the difference of the two reactive distances, i.e., the distance between the oxygen of the nucleophilic water and the γ-phosphorus, d(Pγ - OW), and the scissile bond length, d(Pγ - O). The other reaction coordinate describes the proton transfer from the nucleophilic water to the γ phosphate, mediated by Ser23618,31, with a linear combination of several antisymmetric stretch coordinates18. In the adiabatic mapping calculations, the threshold for the GRMS is set to 0.05 kcal/(mol·Å); test calculation shows that decreasing GRMS to 0.01 kcal/(mol·Å) leads to very similar results. The adiabatic mapping is carried out using B3LYP/6-31G(d)/MM. To refine energetics, single point calculations with a larger basis set (6-311+G*)45,48 are carried out. For comparison, adiabatic mapping calculations have also been carried out at the SCC-DFTBPR/MM level, which gives qualitatively similar results (see Supporting Information).

We note that the mechanism involving Glu459 as the general base7,30 is not addressed here. This is motivated by the observation that Glu459 is engaged in a salt-bridge interaction with Arg238 in the pre-powerstroke state (Fig. 1) and therefore unlikely to be a capable catalytic base. In a separate study (Yang and Cui, manuscript in preparation), however, we have found preliminary evidence supporting the feasibility of this somewhat unexpected mechanism in, at least, the S236A mutant.

Figure 1
Comparison of the optimized active site of myosin motor domain with ATP bound to the (a) open (post-rigor) and (b) closed (pre-powerstroke) conformation state. The geometries have been optimized at the B3LYP/6-31G(d)/MM level (see text). Note the significant ...


A. Anomeric effects and ATP “activation” in myosin

1. ATP geometry in different myosin motor domain conformations

In the pre-powerstroke state, the active site adopts the closed configuration (see Fig. 1), where the most striking feature in the optimized structure is the substantially elongated Pγ-O bond distance in ATP. The value of 1.77 Å is significantly longer than the values of 1.69-1.70 Å found for unprotonated methyl triphosphate in solution by several molecular simulation studies56,57. In fact, this particular feature has been observed in several previous analyses of phosphoryl transfer enzymes and ATPases. The Pγ-O distance in the same myosin-ATP system was optimized to be 1.78 Å in the work of Grigorenko et al.7, who used a somewhat lower level of theory (Hartree-Fock with effective-core-potential and the corresponding double-zeta plus polarization quality bases) in a QM/MM framework. Similarly, ATP with an elongated Pγ-O bond (1.77 Å) was observed in the QM/MM study of F1-ATPase by Dittrich et al.8, although a smaller basis set (6-31G) was used. Finally, calculations of the Ras-GAP-GTP complex6,9 found that the GTP ligand has a long Pγ-O distance of 1.78 Å; interestingly, the distance is slightly shorter (~1.73 Å) in the absence of GAP binding, which seems to be consistent with the role of GAP in activating GTP hydrolysis. Indeed, infrared (IR) spectroscopy studies also supported the notion that the scissile bond becomes weaker upon GTP binding to Ras, especially in the Ras-GAP complex2,3,10.

The somewhat similar finding of the current study is that the Pγ-O distance in ATP is substantially shorter when the myosin active site adopts the open configuration. In the post-rigor state, the QM/MM optimized value for the Pγ-O bond length is 1.72 Å (Fig. 1), which is 0.05 Å shorter than that in the pre-powerstroke state. Therefore, the Pγ-O bond is significantly elongated only when ATP is bound to the hydrolysis-competent state of myosin, which begs the question whether this suggests that ATP is chemically activated and poised to hydrolyze only in the pre-powerstroke state of myosin, similar to the proposed role of GAP in activating GTP in Ras9,10. We note that this structural feature is not directly available from crystallographic studies because the X-ray structure of the pre-powerstroke state could only be captured with a transition state analogue, ADPVO4, as the ligand34.

2. NBO analysis of electronic structure

A well-defined framework for revealing the connection between reactivity and structure is the NBO analysis58,59. Following this framework, Ruben et al.60,61 has demonstrated that the selective destabilization of the scissile bond in “high energy” phosphate compounds is largely due to the generalized anomeric effect, which in electronic structure terms (see Fig. 2) is related to the charge transfer interaction from the lone pair orbital of the non-bridging oxygens (Onon-bridge) to the anti-bonding orbital of the neighboring P-O bond that involves the bridging oxygen (Obridge). Using 10 model phosphoryl compounds as examples, they found strong correlations between experimental hydrolysis free energy, optimized Pγ-O bond distance and second-order perturbation estimates for the interaction between Onon-bridge (O,O and O) lone pairs and the anti-bonding orbital of Pγ-O.

Figure 2
Illustrations of the anomeric effects in tri-phosphate. (a) A schematic illustration for the electron transfer from a lone-pair on the non-bridging oxygen to the anti-bonding orbital of the scissile P-O bond; (b) 3D (left) and 2D (right) representations ...

Encouraged by these findings and using B3LYP/6-311+G* single point calculations at the QM/MM optimized structures, we have carried out NBO analysis to characterize the magnitude of anomeric effects in ATP bound to the open (post-rigor) and closed (pre-powerstroke) configurations of the myosin motor domain (see Table 1). For the closed configuration, second-order perturbation theory analysis of the Fock matrix in the NBO basis finds the interaction energy between the Onon-bridge lone pairs and the anti-bonding orbital of Pγ-O as 21.3, 23.0 and 25.3 kcal/mol, respectively, which are consistent with Ruben et al.'s observation60,61 of ~25 kcal/mol for each such interaction in model phosphates in the gas phase. Those authors also generated a correlation between the magnitude of the charge transfer interactions and the Pγ-O distance; a bond distance of 1.77 Å corresponds to ~72 kcal/mol of charge transfer interactions, which is again close to the sum of interactions found here (~69.6 kcal/mol). Consistent with these significant orbital interactions, the occupation number of the Pγ-O anti-bonding orbital is as high as 0.27e, significantly larger than the values for the Pγ-Onon-bridge anti-bonding orbitals, which are in the range of 0.12-0.14e; these values are similar to the results of NBO analysis for Ras-GAP-GTP9. Similarly, the Wiberg bond order62 for the Pγ-O bond is 0.46, indicating a significantly weakened covalent bond.

Table 1
NBO analysis results for anomeric effects associated with the Pγ-O bond of ATP in the open and closed conformations of the myosin motor domain and of Methyl-triphosphate (MTP) in solutiona

For the open active site configuration (see Table 1), the anomeric interactions become weaker, in line with the shorter Pγ-O distance. Second-order perturbation theory analysis finds the anomeric interactions to be 19.7, 20.6 and 21.8 kcal/mol, respectively, thus the sum of interactions is weaker than that in the closed configuration by approximately 7 kcal/mol; this decrease in the magnitude of anomeric interactions as the Pγ-O distances gets shortened by 0.05 Å is again in nearly quantitative agreement with the correlation found in the study of Ruben and co-workers60. Along the same line, the occupation number of the Pγ-O antibonding orbital is 0.24e and the Wiberg bond order is 0.52.

To identify the origin for the difference in anomeric effects as the active site undergoes the relatively subtle open/close transition (see Fig. 1), perturbation analysis is carried out in which the NBO analysis is repeated after the MM contribution from a specific residue is turned off. The results indicate that the residual contributions are relatively small and do not differ dramatically in the two myosin conformations, which highlights the intramolecular nature of the anomeric effect; e.g., turning off the charges on Arg238 has the largest differential effect, which perturbs the total anomeric effect by 0.15 and 2.18 kcal/mol in the closed and open states, respectively. Therefore, the majority of the difference in the anomeric effect between the open/closed active sites comes from the Pγ-O bond distance variation; this is confirmed by a partial optimization calculation in which the Pγ-O bond in the closed state is restrained to be 1.72 Å, and subsequent NBO analysis finds that the anomeric effect is reduced by 4.1 kcal/mol compared to the unconstrained result. This highlights that it is not straightforward to determine the causal relationship between the magnitude of the anomeric effect and the geometry of the ATP molecule, despite the good correlation between them found here and in previous studies60,61.

3. Energetics of the Pγ-O bond stretch

The geometrical and NBO analysis discussed above clearly show that the ATP in the closed (pre-powerstroke) active site is substantially different from the nucleotide in solution and the open (post-rigor) motor domain. To state that these structural and electronic changes indicate that ATP is significantly “activated” in the closed active site, however, requires a careful characterization of energetics. Accordingly, we have carried out potential energy scan along the Pγ-O bond for an isolated tri-phosphate; for each Pγ-O distance, the rest of structural parameters are optimized in the gas phase, and single point energy calculations are then carried out using a continuum solvent model (PCM). As shown in Fig. 3, the energy changes very modestly as the Pγ-O bond length increases from the reported solution value (~1.69 Å)56,57 through the optimal value in the open motor domain (1.72 Å) to the optimized value in the closed active site (1.77 Å). The energy consequence for lengthening the Pγ-O distance from 1.69 to 1.77 Å is only about 1 kcal/mol and this result is not sensitive to the atomic radii in the PCM model. Therefore, we have to conclude that despite the significant increase in the Pγ-O distance upon binding to the catalytically competent form of myosin, the ATP is not significantly activated compared to solution. Therefore, the dramatic decrease in the hydrolysis barrier in the motor domain14,15,27 has to come largely from the stabilization of the transition state rather than destabilization of the ground state relative to solution.

Figure 3
PCM single point energy (B3LYP/6-311+G(d)) as a function of the Pγ-O bond distance in a methyl-triphosphate molecule; all other degrees of freedom have been optimized in the gas-phase at the B3LYP/6-31G(d) level. The result is not sensitive ...

B. The nature of ATP hydrolysis transition state in myosin

Based on the two-dimensional adiabatic map at the B3LYP/MM level (Fig. 4), a single barrier of ~26 kcal/mol separates the reactant and the intermediate state. Similar to observations from our recent study using SCC-DFTBPR/MM18, the Pγ-O bond is largely broken in the intermediate, with a long distance of 2.84 Å. The approximate saddle point from the adiabatic mapping calculations can be further optimized with the conjugate peak refinement (CPR) method63, which leads to a true saddle point that differs in energy by only ~0.5 kcal/mol; the two reaction coordinate values at the optimized saddle point are -0.47 and 0.12, which can be compared to the approximate saddle point on the 2-dimensional adiabatic map at ~(-0.5,0.2) (see Fig. 4). Therefore, by both geometry and energy criteria, the reaction coordinates chosen for the adiabatic mapping calculations are effective.

Figure 4
Two-dimensional adiabatic map for ATP hydrolysis in myosin at the B3LYP/6-31G*/CHARMM level (with single point B3LYP/6-311+G*/MM energies); the x-axis is defined as the difference between two reactive P-O distances, d(P γ - OW ) - d(P γ ...

We note that the potential energy barrier of ~ 26 kcal/mol is substantially higher than the free energy barrier estimated based on experimental rate constant14,27, which is ~15-17 kcal/mol. This discrepancy, however, is rather well understood based on our recent SCCDFTBPR/MM studies. As discussed in Ref.18, SCC-DFTBPR/MM minimum energy path (MEP) calculations also found a rather high barrier similar to the current adiabatic mapping calculations and our previous MEP analysis16 at the B3LYP/6-31G(d,p)//HF/3-21+G level. By contrast, potential of mean force (PMF) simulations at the SCC-DFTBPR/MM level found a free energy barrier very close to the experimental value. Comparison of structures from the MEP and PMF simulations found that the majority of residues/water molecules in the active site have similar structures, except for Ser181 and water molecules coordinated to its sidechain. As the γ phosphate gets protonated during the hydrolysis, it is energetically favorable for Ser181 to rotate and form alternative hydrogen-bonding interactions; since such an isomerization is coupled to reorientation of nearby water molecules, it was sampled during PMF simulations but not in the MEP calculations. Perturbative analysis confirmed18that differences in interactions due to Ser181 and nearby water account for most discrepancy between the MEP and PMF barriers. Since other structural features of the active site are highly consistent between the PMF and MEP simulations, it is meaningful to use the MEP (and adiabatic mapping) results to further explore the properties of the transition state region.

The most striking feature of the two-dimensional adiabatic map is that the saddle point region is very flat energetically (Fig.4). For a rather large region that covers from -1.1Å to 0.3Å in the x-axis and from -0.2Å A to 0.7Å A in the y-axis (indicated by a rectangle in Fig. 4), the average energy variation is only 1.8 kcal/mol. More importantly, when examining the structures in this region, various combinations of Pγ-O and Pγ-OW distances that correspond to different classical mechanisms can be observed (Fig. 5). For example, in structure (0.285, 0.382), the Pγ-O and Pγ-OW distances are 1.86 and 2.15 Å, respectively, which correspond to an associative-like “transition state”. However, in structure (-0.006, 0.025), the Pγ-O and Pγ-OW distances are both 2.01 Å, which show significant concerted feature. In structure (-0.297, 0.382), the Pγ-O and Pγ-OW distances are 2.19 Å and 1.89 Å A respectively, reminiscent of a dissociative “transition state”. Importantly, the energies for the three structures are 25.3, 23.8 and 25.8 kcal/mol, respectively, and therefore are nearly degenerate. To further confirm that the observed small energy differences are not an artifact due to the neglect of Ser181 reorientation in the current calculations, we examine both NBO charges on ATP (see Table 2) and the energetic contribution from Ser181 in the four structures shown in Fig. 5. As shown in Table 2, the only significant trend is that the charge on O becomes increasingly more negative as the structure becomes more dissociative in nature, although the magnitude of variation is rather small (~0.05 e). Therefore, it is unlikely that interaction with the enzyme environment has a very significant impact on the relative stability of the four representative “transition-state-like” structures. Indeed, the perturbative energetic contribution from Ser181 in the four structures, relative to the reactant state, is -3.1, -5.2, -3.8 and -3.8 kcal/mol, respectively. The slightly higher unfavorable contribution from Ser181 in the “associative” structure (Structure 2 in Fig. 5) is qualitatively consistent with the observation that the PMF saddle point at the SCC-DFTBPR/MM level is more associative-like (i.e., a shorter Pγ-O bond) than the MEP saddle points (see Table 2 in the Supporting Information in Ref.18).

Figure 5
Structures near the saddle point (shown in a) that are representative of associative, concerted and dissociative mechanisms, respectively; they have the coordinates of (b) (0.285, 0.382) (c) (-0.006, 0.025) and (d) (-0.297, 0.382) on the two-dimensional ...
Table 2
Natural Population Chargesa on key atoms in ATP in the four structures shown in Fig. 5.b

In short, the two-dimensional adiabatic map results raise the issue that the energetic difference between reaction mechanisms following the classical nomenclature, at least for ATP hydrolysis in myosin, can be rather small. Various structures found in the saddle point region, which could be classified to represent different mechanisms, have very similar energies and therefore all thermally accessible; this needs to be further corroborated using more extensive explorations of the free energy surface beyond a few stationary points. The lack of a clear distinction between different extreme types of transition states as envisioned in classical nomenclature22,23 has been discussed for the hydrolysis of phosphate compounds in the aqueous solution24,26,64, and our study highlights this point again for a typical ATPase.


In the analysis of factors that contribute to enzyme catalysis, it is generally useful to employ a multi-facted approach in which multiple properties of the substrate and substrate-enviornment interactions are analyzed. From this perspective, spectroscopic analysis such as infrared and NMR spectroscopies can potentially reveal structural differences when the substrate is in different environments. Similarly, electronic structure techniques such as NBO analysis can quantify changes in the electronic structure of the substrate as its environment varies, and therefore providing more concrete observables for the discussion of mechanistic hypothesis such as “ground state destabilization”. On the other hand, it is useful to note that the potential caveat is that the observables in these analyses tend to very sensitive to perturbations in the substrate. In other words, notable shifts in spectra or NBO results may not necessarily correlate to significant energetic perturbations. Take the system in the current study as an example, a stretch of ~0.08 Å in the covalent Pγ-O bond and 7 kcal/mol in the anomeric effect are rather significant changes in structural and NBO terms. However, these changes are correlated to an increase in energy of about 1 kcal/mol in the ATP molecule, which is fairly insignificant in the context of enzyme catalysis. Therefore, we have to conclude that ground state destabilization is not a significant contribution to the hydrolysis of ATP in myosin, and it is likely that this conclusion also applies to GTP hydrolysis in Ras-GAP, for which previous studies focused on the charge distribution9,10 without an explicit analysis of energetics associated with Pγ-O elongation. Our study emphasizes that spectroscopic, structural and electronic structural analyses, as valuable as they are, should always be complemented with an evaluation of the relevant energetic contribution for the purpose of determining whether a factor is essential to catalysis.

For defining the nature of the transition state and therefore the reaction mechanism, it is typical to focus on the properties of a single structure, such as the saddle point optimized on the potential energy surface or, in some cases, on the free energy surface65. This is meaningful only when representative structures from different mechanisms have distinct (free) energies. The current study clearly indicates that this may not be the case in some systems. If we only examine the structural features of the saddle point, the hydrolysis of ATP appears to be associative in nature. By exploring the potential energy surface around the saddle point, however, we realize that structures that have clear dissociative and concerted features are in fact very close (within 0.5-1 kcal/mol) in energy. Therefore, one can't conclude at all that other pathways are much less favorable than the associative mechanism. Whether this translates to other phosphate hydrolysis systems remains to be carefully analyzed. However, our observation highlights the importance of looking beyond the saddle point(s) (even on free energy surface) for characterizing the nature of transition state and reaction mechanism, which is likely valid in general. Along this line, broader explorations of the energy landscape that involve either mapping multi-dimensional (free) energy surfaces66 or transition state ensemble analysis using transition state sampling67,68 are required. The high computational cost associated with these types of calculations underlines the value of developing reliable semi-empirical QM methods69-73 for an in-depth mechanistic analysis of complex systems.

Supplementary Material



The research discussed here has been supported from the National Institutes of Health (R01-GM071428). Computational resources from the National Center for Supercomputing Applications at the University of Illinois are greatly appreciated.


Supporting Information Available Adiabatic Mapping calculations at the SCC-DFTBPR/MM level are compared to the B3LYP/MM calculations. This material is available free of charge via the Internet at


1. Alberts B, Bray D, Lewis J, Raff M, Roberts K, Watson JD. Molecular biology of the cell. Garland Publishing, Inc.; 1994.
2. Allin C, Ahmadian MR, Wittinghofer A, Gerwert K. Monitoring the gap catalyzed HRas GTPase reaction at atomic resolution in real time. Proc. Acad. Natl. Sci. USA. 2001;98:7754–7759. [PubMed]
3. Klahn M, Schlitter J, Gerwert K. Theoretical IR spectroscopy based on QM/MM calculations provides changes in charge distribution, bond lengths, and bond angles of the GTP ligand induced by the Ras-protein. Biophy. J. 2005;88:3829–3844. [PubMed]
4. Barth A. Infrared spectroscopy of proteins. Biochim. Biophys. Acta - Bioeng. 2007;1767:1073–1101. [PubMed]
5. Barth A, Bezlyepkina N. P-O bond destabilization accelerates phosphoenzyme hydrolysis of sarcoplasmic reticulum Ca2+-ATPase. J. Biol. Chem. 2004;179:51888–51896. [PubMed]
6. Grigorenko BL, Nemukhin AV, Topol IA, Cachau RE, Burt SK. QM/MM modeling the Ras-GAP catalyzed hydrolysis of guanosine triphosphate. Proteins: Struct., Funct., and Bioinf. 2005;60:495–503. [PubMed]
7. Grigorenko BL, Rogov AV, Topol IA, Burt SK, Martinez HM, Nemukhin AV. Mechanism of the myosin catalyzed hydrolysis of atp as rationalized by molecular modeling. Proc. Acad. Natl. Sci. U.S.A. 2007;104:7057–7061. [PubMed]
8. Dittrich M, Hayashi S, Schulten K. ATP hydrolysis in the β(TP) and β(DP) catalytic sites of F1-ATPase. Biophy. J. 2004;87:2954–2967. [PubMed]
9. Grigorenko BL, Nemukhin AV, Shadrina MS, Topol IA, Burt SK. Mechanisms of Guanosine Triphosphate hydrolysis by Ras and Ras-GAP proteins as rationalized by ab initio QM/MM simulations. Proteins: Struct., Funct., and Bioinf. 2007;66:456–466. [PubMed]
10. te Heesen H, Gerwert K, Schlitter J. Role of the the arginine finger in Ras·RasGAP revealed by QM/MM calculations. FEBS Lett. 2007;581:5677–5684. [PubMed]
11. Vale RD, Milligan RA. The way things move: looking under the hood of molecular motor proteins. Science. 2000;288:88–95. [PubMed]
12. Eisenberg E, Hill TL. Muscle contraction and free energy transduction in biological systems. Science. 1985;227:999–1006. [PubMed]
13. Hill TL. Free energy transduction in biology. Academic Press; New York: 1977.
14. Geeves MA, Holmes KC. Structural mechanism of muscle contraction. Annu. Rev. Biochem. 1999;68:687–728. [PubMed]
15. Geeves MA, Holmes KC. The molecular mechanism of muscle contraction. Adv. Prot. Chem. 2005;71:161–193. [PubMed]
16. Li GH, Cui Q. Mechanochemical coupling in myosin: A theoretical analysis with molecular dynamics and combined QM/MM reaction path calculations. J. Phys. Chem. B. 2004;108:3342–3357.
17. Yu H, Ma L, Yang Y, Cui Q. Mechanochemical coupling in myosin motor domain, I. Equilibrium active site simulations. PLoS Comput. Biol. 2007;3:0199. [PMC free article] [PubMed]
18. Yang Y, Yu H, Cui Q. Extensive conformational changes are required to turn on atp hydrolysis in myosin. J. Mol. Biol. 2008;381:1407–1420. [PMC free article] [PubMed]
19. Cleland WW, Hengge AC. Enzymatic mechanisms of phosphate and sulfate transfer. Chem. Rev. 2006;106:3252–3278. [PubMed]
20. Warshel A, Sharma PK, Kato M, Xiang Y, Liu HB, Olsson MHM. Electrostatic basis for enzyme catalysis. Chem. Rev. 2006;106:3210–3235. [PubMed]
21. Yang Y, Cui Q. Interactions between phosphate and water in solution: A Natural Bond Orbital based analysis in a QM/MM framework. J. Phys. Chem. B. 2007;111:3999–4002. [PMC free article] [PubMed]
22. Guthrie RD, Jencks WP. IUPAC recommendations for the representation of reaction mechanisms. Acc. Chem. Res. 1989;22:343–349.
23. Hassett A, Blättler W, Knowles JR. Pyruvate kinase: Is the mechanism of phospho transfer associative or dissociative? Biochem. 1982;21:6335–6340. [PubMed]
24. Kamerlin SCL, Florian J, Warshel A. Associative versus dissociative mechanisms of phosphate monoester hydrolysis: On the interpretation of activation entropies. ChemPhysChem. 2008;9:1767–1773. [PubMed]
25. Admiraal S, Herschlag D. Mapping the transition state for atp hydrolysis: implications for enzymatic catalysis. Chem. Biol. 1995;2:729–739. [PubMed]
26. Aqvist J, Kolmodin K, Florian J, Warshel A. Mechanistic alternatives phosphate monoester hydrolysis: what conclusions can be drawn from available experimental data? Chem. & Biol. 1999;6:R71–R80. [PubMed]
27. Trentham DR, Eccleston JF, Bagshaw CR. Kinetic-analysis of ATPase mechanisms. Q. Rev. Biophys. 1976;9:217–281. [PubMed]
28. Deng H, Wang JH, Callender RH, Grammer JC, Yount RG. Raman difference spectroscopic studies of the myosin S1·MgADP vanadate complex. Biochem. 1998;37:10972–10979. [PubMed]
29. Schwarzl SM, Smith JC, Fischer S. Insights into the chemomechanical coupling of the myosin motor from simulation of its ATP hydrolysis mechanism. Biochem. 2006;45:5830–5847. [PubMed]
30. Onishi H, Ohki T, Mochizuki N, Morales MF. Early stages of energy transduction by myosin: Roles of Arg in switch I, of Glu in switch II, and of the salt-bridge between them. Proc. Natl. Acad. Sci. U.S.A. 2002;99:15339–15344. [PubMed]
31. Rayment I. The structural basis of the myosin ATPase activity. J. Biol. Chem. 1996;271:15850–15853. [PubMed]
32. Gulick AM, Bauer CB, Thoden JB, Rayment I. X-ray structure of the MgADP, MgATPγS, and MgAMPPNP complexes of the dictostelium discoideum myosin motor domain. Biochem. 1997;36:11619–11628. [PubMed]
33. Bauer CB, Holden HM, Thoden JB, Smith R, Rayment I. X-ray structures of the Apo and MgATP-bound states of Dictyostelium discoideum myosin motor domain. J. Biol. Chem. 2000;275:38494–38499. [PubMed]
34. Smith CA, Rayment I. X-ray structure of the Magnesium(II) · ADP · Vanadate complex of the Dictyostelium discoideum myosin motor domain to 1.9 Å resolution. Biochemistry. 1996;35:5404–5417. [PubMed]
35. Warshel A. Computer Modeling of Chemical Reactions in Enzymes and Solution. Wiley; New York: 1991.
36. Field MJ, Bash PA, Karplus M. A combined quantum-mechanical and molecular mechanical potential for molecular-dynamics simulations. Journal of Computational Chemistry. 1990;11(6):700–733.
37. Gao J. In: Reviews in Computational Chemistry VII. Lipkowitz KB, Boyd DB, editors. VCH; New York: 1995.
38. A. D. M., Bashford D, Bellott M, R. L. D., Evenseck JD, Field MJ, Fischer S, Gao J, Guo H, Ha S, Joseph-McCarthy D, Kuchnir L, Kuczera K, Lau FTK, Mattos C, Michnick S, Ngo T, Nguyen DT, Prodhom B, W. E. R., Roux B, Schlenkrich M, Smith JC, Stote R, Straub J, Watanabe M, Wiórkiewicz-Kuczera J, Yin D, Karplus M. All-atom empirical potential for molecular modeling and dynamics studies of proteins. Journal of Physical Chemistry B. 1998;102:3586–3616. [PubMed]
39. Elstner M, Porezag D, Jungnickel G, Elsner J, Haugk M, Frauenheim T, Suhai S, Seifert G. Self-Consistent-Charge Density-Functional Tight-Binding method for simulations of complex materials properties. Phys. Rev. B. 1998;58(11):7260–7268.
40. Yang Y, Yu H, York D, Elstner M, Cui Q. Description of phosphate hydrolysis reactions with the Self-Consistent-Charge Density-Functional-Tight-Binding (SCC-DFTB) theory 1. Parametrization. J. Chem. Theo. Comp. 2008;4:2067–2084. [PMC free article] [PubMed]
41. Becke AD. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A. 1988;38:3098–3100. [PubMed]
42. Becke AD. Density-functional thermochemistry. iii. the role of exact exchange. J. Chem. Phys. 1993;98:5648–5652.
43. Lee C, Yang W, Parr RG. Development of the colle-salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B. 1988;37:785–789. [PubMed]
44. Haharan PC, Pople JA. The influence of polarization functions on molecular orbital hydrogenation energies. Theor. Chim. Acta. 1973;28:213–222.
45. Krishnan R, Binkley JS, Seeger R, Pople JA. Self-consistent molecular orbital methods. xx. a basis set for correlated wavefunctions. j. Chem. Phys. 1980;72:650–654.
46. Frisch MJ, Pople JA, Binkley JS. Self-consistent molecular-orbital methods .25. supplementary functions for gaussian-basis sets. J. Chem. Phys. 1984;80:3265–3269.
47. Clark T, Chandrasekhar J, Spitznagel GW, Schleyer P. v. R. Efficient diffuse function-augmented basis-sets for anion calculations .3. the 3-21+g basis set for 1st-row elements, li-f. J. Comp. Chem. 1983;4:294–301.
48. McLean AD, Chandler GS. Contracted gaussian-basis sets for molecular calculations .1. 2nd row atoms, Z=11-18. J. Chem. Phys. 1980;72:5639–5648.
49. Brooks BR, Bruccoleri RE, Olafson BD, States DJ, Swaminathan S, Karplus M. Charmm: A program for macromolecular energy, minimization and dynamics calculations. Journal of Computational Chemistry. 1983;4(2):187–217.
50. Shao Y, Fusti-Molnar L, Jung Y, Kussmann J, Ochsenfeld C, Brown ST, Gilbert ATB, Slipchenko LV, Levchenko SV, O'Neill DP, Distasio RA, Jr., Lochan RC, Wang T, Beran GJO, Besley NA, Herbert JM, Lin CY, Van Voorhis T, Chien SH, Sodt A, Steele RP, Rassolov VA, Maslen PE, Korambath PP, Adamson RD, Austin B, Baker J, Byrd EFC, Dachsel H, Doerksen RJ, Dreuw A, Dunietz BD, Dutoi AD, Furlani TR, Gwaltney SR, Heyden A, Hirata S, Hsu C-P, Kedziora G, Khalliulin RZ, Klunzinger P, Lee AM, Lee MS, Liang W, Lotan I, Nair N, Peters B, Proynov EI, Pieniazek PA, Rhee YM, Ritchie J, Rosta E, Sherrill CD, Simmonett AC, Subotnik JE, Woodcock HL, III, Zhang W, Bell AT, Chakraborty AK, Chipman DM, Keil FJ, Warshel A, Hehre WJ, Schaefer HF, III, Kong J, Krylov AI, Gill PMW, Head-Gordon M. Advances in methods and algorithms in a modern quantum chemistry program package. Phys. Chem. Chem. Phys. 2006;8:3172–3191. [PubMed]
51. Woodcock L, Hodoek M, Gilbert ATB, Gill PMW, Schaefer HF, Brooks BR. Interfacing q-chem and charmm to perform qm/mm reaction path calculations. J. Comp. Chem. 2007;28:1485–1502. [PubMed]
52. Glendening ED, Badenhoop JK, Reed AE, Carpenter JE, Bohmann JA, Morales CM, Weinhold F. NBO 5.0. 2001. nbo5.
53. Cossi M, Barone V, Cammi R, Tomasi J. Ab initio study of solvated molecules: A new implementation of the polarizable continuum model. Chem. Phys. Lett. 1996;255:327–335.
54. Barone V, Cossi M, Tomasi J. A new definition of cavities for the computation of solvation free energies by the polarizable continuum model. J. Chem. Phys. 1997;107:3210–3221.
55. Brooks CL, III, Karplus M, Pettitt BM. Proteins: A theoretical perspective of dynamics, structure, and thermodynamics. Wiley and Sons; New York: 1988.
56. Grigorenko BL, Rogov AV, Nemukhin AV. Mechanism of triphosphate hydrolysis in aqueous solution: QM/MM simulations in water clusters. J. Phys. Chem. B. 2006;110:4407–4412. [PubMed]
57. Akola J, Jones RO. Atp hydrolysis in water - a density functional study. J. Phys. Chem. B. 2003;107:11774–11783.
58. Reed AE, Curtiss LA, Weinhold F. Intermolecular interactions from a natural bond orbital, donor-acceptor viewpoint. Chem. Rev. 1988;88:899–926.
59. Weinhold F, Landis CR. Valency and Bonding. Cambridge University Press; 2005.
60. Ruben EA, Plumley JA, Chapman MS, Evanseck JD. Anomeric effect in “high energy” phosphate bonds. selective destabilization of the scissile bond and modulation of the exothermicity of hydrolysis. J. Am. Chem. Soc. 2008;130:3349–3358. [PubMed]
61. Ruben EA, Chapman MS, Evanseck JD. Generalized anomeric interpretation of the “high-energy” N-P bond in N-methyl-N'-phosphorylguanidine: Importance of reinforcing stereoelectronic effects in “high-energy” phosphoester bonds. J. Am. Chem. Soc. 2005;127:17789–17798. [PubMed]
62. Wiberg KB. Application of pople-santry-segal cndo method to cyclopropylcarbinyl and cyclobutyl cation and to bicyclobutane. Tetrahedron. 1968;24:1083.
63. Fischer S, Karplus M. Conjugate peak refinement: an algorithm for finding reaction paths and accurate transition states in systems with many degrees of freedom. Chem. Phys. Lett. 1992;194:252–261.
64. Florian J, Warshel A. Phosphate ester hydrolysis in aqueous solution: Associative versus dissociative mechanisms. J. Phys. Chem. B. 1998;102:719–734.
65. Hu H, Lu ZY, Yang WT. Qm/mm minimum free-energy path: Methodology and application to triosephosphate isomerase. J. Chem. Theo. Comp. 2007;3:390–406. [PMC free article] [PubMed]
66. Barducci A, Bussi G, Parrinello M. Well-tempered metadynamics: A smoothly converging and tunable free-energy method. Phys. Rev. Lett. 2008;100 020603. [PubMed]
67. Bolhuis PG, Chandler D, Dellago C, Geissler PL. Transition path sampling: Throwing ropes over rough mountain passes, in the dark. Annu. Rev. Phys. Chem. 2002;53:291–318. [PubMed]
68. Crehuet R, Field MJ. A transition path sampling study of the reaction catalyzed by the enzyme chorismate mutase. J. Phys. Chem B. 2007;111:5708–5718. [PubMed]
69. Thiel W. Perspectives on semiempirical molecular orbital theory. Adv. Chem. Phys. 1996;93:703–757.
70. Repasky MP, Chandrasekhar J, Jorgensen WL. PDDG/PM3 and PDDG/MNDO: Improved semiempirical methods. J. Comp. Chem. 2002;23:1601–1622. [PubMed]
71. Nam K, Cui Q, Gao J, York DM. A specific reaction parameterization for the am1/d hamiltonian for transphosphorylation reactions. J. Theo. Comp. Chem. 2007;3:486–504.
72. Elstner M. The scc-dftb method and its application to biological systems. Theo. Chem. Acc. 2006;116:316–325.
73. Yang Y, Yu H, York D, Cui Q, Elstner M. Extension of the self-consistent-charge tightbinding-density-functional (scc-dftb) method: third order expansion of the dft total energy and introduction of a modified effective coulomb interaction. J. Phys. Chem. A. 2007;111:10861–10873. [PubMed]