In the QM/MM calculations, the QM system was defined as the phosphate group () and the attacking water molecule. shows four distances relevant to the phosphate hydrolysis mechanism. The initial reaction coordinate Qe= r1− r2 was chosen as the difference in length of the breaking and forming bonds, i.e., the phosphorus atom and the leaving O3' oxygen atom (r1), and the oxygen of the attacking water molecule and the phosphorus atom (r2). Qe probes primarily the electron-transfer part of the hydrolysis reaction. This choice of reaction coordinate allows us to distinguish between dissociative and associative paths. As the reaction is forced to progress by advancing the windows in the US simulations, the two extreme responses of the system are either an increase in the distance between the phosphorus and the O3' oxygen of the leaving group (which would lead to a dissociative mechanism), or a decrease in the distance between the attacking water oxygen and the phosphorus.
To initialize the US/PMF simulations, we prepared starting structures of the system through a series of energy minimizations along the reaction coordinate Qe. Our preliminary QM/MM minimization studies showed hysteresis when the Qe reaction coordinate was used (data not shown). We started minimizations from the reactant state and moved towards the product state and vice versa. We found that the forward and backward minimum energy pathways differed at the transition state region, although the endpoints (RS and PS) were identical.
The PMF profile obtained from the standard US simulations along Qe
appears well converged, as shown in . The individual free energy curves for different US windows overlap, have similar local curvatures, and the overall free energy profile appears well connected. The PMF suggests an associative mechanism with a low barrier for the first step of the hydrolysis reaction, resulting in a stable penta-coordinated intermediate. The barrier for this step is about 9.5 kcal/mol, significantly lower than the experimental estimate of about 21 kcal/mol35
Figure 3 One-dimensional free energy profiles of RNase H reaction. (a) Free energy profile along the reaction coordinate Qe obtained from PMF/US simulations. The overall free energy profile (circles) was obtained by using WHAM43. The lines show the unbiased free (more ...)
Selected average atomic distances, r1, r2, r3 and r4 () observed along Qe are shown in (panels a and c for standard US and RE/US, respectively). Consistent with an associative mechanism, the distance between the water molecule and the phosphate (r2) decreases, while the leaving-group distance (r1) remains nearly constant as the reaction progresses. However, other relevant distances, such as r3 and r4 (involving the hydrogen atom of the attacking water molecule) present discontinuous jumps along Qe.
Average atomic distances r1, r2, r3 and r4 () observed along the reaction coordinates Qe (panels a and c) and Qep (panels b and d). Results from (a,b) US and (c,d) RE/US.
Such discontinuities in distances not included in the reaction coordinate Qe can artificially lower the apparent free energy barrier. Thus, we reexamined the free energy simulations at the transition state region for discontinuities, following the arguments given in the discussion of . To systematically search for additional important coordinates, we developed the following simple analysis. We calculated the averages and the standard deviations of the atom-atom distances that involved each of the QM atoms and all atoms within 4.0 Å distance of the QM subsystem (a total of about 1000 atom pair distances, dij). For each of the neighboring US windows, we analyzed the difference between the average dij values and considered the path to be continuous if the difference did not exceed twice the sum of the standard deviations of dij in the two neighboring windows. By using this procedure, we discovered that several distances at the transition state were not continuous. The majority of these distances involved the hydrogen atom of the attacking water molecule that was transferred to the O1P oxygen of the phosphate group in a sharp, jump-like transition (curves for r3 and r4 in , with r3 being the distance between the hydrogen and the oxygen of the attacking water, and r4 the distance between the same hydrogen and the O1P oxygen atom of the phosphate).
To include these missing degrees of freedom, we extended the reaction coordinate Qe by combining it with Qp = r3 − r4 in a new coordinate Qep = Qe + Qp. The combined coordinate Qep accounts for the forming and breaking bonds by adding the distance of the hydrogen and the oxygen of the water molecule (r3) and by subtracting the distance between the hydrogen of the water and the oxygen of the phosphate (r4) ().
We performed a new set of US/PMF simulations along Qep with 33 windows along this reaction coordinate. The centers of the harmonic biasing potentials were equally spaced from Qep = −2.6 Å to 0.6 Å, with a force constant of 2000 kcal/(molÅ2). In addition to the umbrella potentials on the new coordinate Qep, we also applied a harmonic potential at a constant value of Qe=−0.70 Å, with a 100 kcal/(molÅ2) force constant to keep the geometries near the transition state along Qe. RE/US simulations were also carried out. We used a force constant of 100 kcal/(molÅ2) for the biasing potential centered at Qep between −1.75 Å and 0.15 Å in 0.05 Å increments. In the US simulations biasing Qep, we found that all distances varied continuously, including r3 and r4, as shown in .
Data from all four simulations were combined to determine the two-dimensional free energy surface in the Qe
plane. As shown in , our best estimate of the barrier of the first step of the hydrolysis reaction is about 22 kcal/mol, in excellent agreement with the experimental value of 20.5 kcal/mol35
, assuming a pre-factor of k0
=1/ps. Importantly, we found that extending the 1D reaction coordinate from Qe
increased the barrier by more than 10 kcal/mol. We note that an earlier study30
also had used Qe
, and found a barrier that was too low compared to experiment.
Figure 5 Two-dimensional free energy surface for the first step of the hydrolysis reaction of RNase H. All simulation data (both standard US and RE/US, with either Qe or Qep as the reaction coordinate) were included when determining the free energy surface (kcal/mol) (more ...)
The projected one-dimensional free energy profiles along Qe (red, bottom axis) and Qep (black, top axis) are shown in . Note that the barrier shown in the one-dimensional profile for Qe is consistently lower than that for Qep, and it is in quantitative agreement with the profile in obtained by using US along Qe alone. We also consistently find that the intermediate state formed as the product of the first step of the reaction has a free energy about 1 kcal/mol below that of the reactant state.
Improvements over standard US often come from using RE/US, which provides enhanced sampling by accelerating equilibration without adding extra cost to the calculation. Compared to standard US, the average values of r1 to r4 in the RE/US simulations follow the reaction coordinate more continuously as the bias changes (). Nevertheless, replica exchange alone could not account for the missing degrees of freedom when only Qe was used, as the discontinuity remained between the US windows at the transition state region for the r3 and r4 distances. By sampling the transition state region using the improved reaction coordinate Qep, we can see in that proton transfer (Qep = −0.8 Å) precedes the electron transfer between the hydroxide and the phosphate (Qep = −0.2 Å, at the barrier top).
The results obtained by searching along Qep
show a clear improvement in calculating the free energy profile. However, as shown in , the orientation of the saddle is not parallel to Qep
in the Qe
plane, and Qep
is not the optimal one-dimensional choice. Ignoring anisotropy in the dynamics, an ideal 1D reaction coordinate for a single saddle point is given by the eigenvector of the second-derivative matrix of the free energy, for which the corresponding eigenvalue is negative. In practice, however, it is often not feasible to evaluate multidimensional Hessian matrices. Here, we note that the r1
) and r3
) coordinates can be considered as independent, and we search for a new, optimized reaction coordinate, Qopt
, as their best linear combination by varying the angle θ in
We integrate the 2D free energy surface to obtain the projected one-dimensional free energy profile along Qopt
and determine the corresponding barrier height as a function of the angle θ. In , the Qe
coordinate is represented by θ = 0° and the corresponding barrier is 9.5 kcal/mol, whereas Qep
corresponds to θ = 45°, with a barrier height of 20.6 kcal/mol. The highest barrier (Qopt
) is at an angle of about θ = 20° with a height of 21.6 kcal/mol, our best estimate of a 1D free energy barrier. We note, however, that the inclusion of additional degrees of freedom beyond Qe
can further increase this barrier, even within the assumptions of our QM/MM model. These results show that at the transition state of the catalytic reaction the optimal coordinate consists of about 30% “proton transfer” and 70% “electron transfer”, which allows us to define a quantitative measure of the coupling between proton transfer and electron transfer processes.
Optimal 1D reaction coordinate Qopt (θ) = Qe cosθ+Qp sinθ. Integrated 1D free energy barriers are calculated as a function of θ by numerical integration using the 2D free energy surface of .
One of the limitations of our method to improve the reaction coordinate is that not in all cases will atomic pair distances be sufficient to capture the important degrees of freedom; instead, collective coordinates may need to be considered44
. Although our method can be useful in diagnosing whether a proposed collective coordinate could be important or not, we have no means of a priori list and test all possible collective coordinates. Another limitation is concerned with the sensitivity of our analysis. If the US windows are very closely spaced, one may observe quasi-continuous trajectories for the average values of the coordinates even if the underlying distribution is bistable. In such cases it will be useful to obtain more detailed information about the actual distribution of the proposed reaction coordinate. In particular, a spike in the variance of a proposed additional coordinate along the reaction could indicate bi-stability.
With regards to the RNase H reaction studied here, we caution against over-interpreting the simulation results despite the apparent agreement with experiment. Besides the choice of the QM system and the reaction coordinate, the size of the basis set and the level of theory also affect the barrier height. To extend the QM/MM model, and to give further details of the complete enzymatic reaction, the dissociation of the leaving sugar group has to be simulated as well. With the current setup, we found that a very high barrier is obtained for the second step, unless additional water molecules or protein atoms are included in the quantum region. Thus for accurate mechanistic details more simulations are underway with an increased size of the QM system.