Properties of dianionic oxyphosphorane intermediates from hybrid QM/MM simulation: Implications for ribozyme reactions

Centre for Structural Biochemistry, Karolinska Institutet, NOVUM, S 141 57 Huddinge, Sweden

Keywords: Hybrid QM/MM, molecular dynamics, free energy, umbrella sampling, transition state, ribozyme

Abstract

Dianionic oxyphosphorane intermediates have been suggested as good models for studying the mechanisms of self-cleavage reactions undergone by some RNA molecules (the so called RNA enzymes or ribozymes). We shall report here the results of a hybrid quantum and molecular mechanical (QM/MM) free energy calculation for obtaining the reaction free energy profiles for such a system. The reactions are simulated in the gas phase as well as within a water sphere in the presence and absence of Mg²⁺ ion. The semiempirical quantum mechanical AM1 hamiltonian describes the oxyphosphorane species and the molecular mechanical CHARMM force field describes the solvent and Mg²⁺. We find that when compared with accurate ab initio calculations, the present method shows a tendency to overestimate the stability of the intermediates. On the other hand, it correctly reproduces the qualitative features, namely, the stabilization of the intermediate with solvation and destabilization in the presence of Mg²⁺. It also shows a drastic lowering of the barrier for bond cleavage in the presence of Mg²⁺ in agreement with ab initio and experimental results. We discuss possible means for increasing the accuracy of the method for subsequent application to model the reactions within ribozymes.

Introduction

Since the discovery of the catalytic power of Tetrahymena pre-rRNA for transesterification reactions of phosphodiester bonds, [reviewed in 1] i.e., cleavage and rejoining of its own nucleotide linkages, the list of such RNA enzymes (ribozymes) is growing. Divalent Mg²⁺ ion is commonly found to be indispensable as a cofactor for these reactions. Such enzymes are thought to be relics of the "RNA world" and as such attract considerable interest. However, the possibility of their practical application is also enormous.

In Tetrahymena-type ribozyme reactions, transesterification is initiated by the attack of the 3'-OH group of the bound guanosine on a phospho-diester linkage [2]. The reaction proceeds with inversion of configuration at the phosphorus centre and thus the 'in-line S_N2(P) mechanism' is the most likely one. [3, 4] Consequently, a pentacoordinate oxyphosphorane intermediate/transition state is postulated.

Quite a few ab initio investigations have been carried out employing oxyphosphorane compounds to gain insight about the reaction mechanisms of ribozymes [5-9]. The problem with these works is that the effect of the bulk environment is not taken into account at all. We are approaching the problem from a different perspective. Our ultimate objective is to model the reaction within the ribozyme. Rigorous quantum calculations are not feasible with such a large molecule, although recently an approximate method, where the contribution from the rest of the enzyme other than the active site had been taken care of using an effective fragment potential, was proposed [10]. The method we are using is generally termed in the current literature as hybrid quantum mechanical molecular mechanical (hybrid QM/MM) method. In this method, the system is partitioned into a classical and a quantum region. The energies and the forces on the quantum atoms are obtained from the expectation values of the QM hamiltonian and its derivatives. The quantum hamiltonian includes electrostatic and van der Waal's interactions with classical atoms. Running molecular dynamics with this method, one can calculate ensemble-averaged thermodynamic quantities, such as free energy change in a chemical reaction.

Several hybrid QM/MM models have been advanced in the last few decades. Molecular mechanical force fields have been combined with semiempirical [11-16], density functional [17], valence bond [18,19] or ab initio Hartree-Fock [20] methods. Their application focuses mainly on solvation phenomena [for a review see 21], while investigation of biochemical problems are also becoming popular [22-27].

Employing a very rigorous calculational procedure for the quantum subsytem, although undoubtedly the most accurate method, is not practicable for a molecular dynamics calculation involving a macromolecule. One has to look for a less time consuming alternative. The semi-empirical hamiltonians are suitable from that point of view. However, a problem with them is that they can be quite unreliable for a particular type of reaction. It is necessary to check their behaviour thoroughly before embarking on a full-fledged simulation with the actual enzymic system. In this report, we present such a test of the semi-empirical AM1 hamiltonian [28] with respect to the phosphodiester bond cleavage reaction using a dianionic trimethoxyphosphorane compound as model (figure 1) [7].

Fig.1 The pentacoordinated oxyphosphorane model

Computational Details

The initial optimised geometry of the intermediate was generated using the AM1 hamiltonian in MOPAC6.0 [29] without imposing any symmetry constraint. This geometry (figure 2) was used as the starting structure for further molecular dynamics calculations.

Fig.2 Initial geometry

Hybrid QM/MM runs were performed as implemented in CHARMM [30]. The dianionic oxyphosphorane species was treated quantum mechanically and the surrounding water molecules and Mg²⁺ ion were treated classically. Umbrella sampling procedure [31] was employed to map the free energy change along the reaction coordinate. The ratio of P-O_ax1 distance and the P-O_ax2 distance was chosen as the reaction coordinate for this study. The final free energy profiles were obtained using the weighted histogram analysis method [32] to stitch together the energy profiles obtained from different windows.

Molecular dynamics runs were performed at a number of windows with suitable umbrella potentials applied on the reaction coordinate of each window. For the gas phase calculations, a total of 24 windows were used. The initial structure was heated from 100 K to 300 K for 1 ps and then equilibrated for 1 ps at 300 K with a value of the reaction coordinate set to 1. The resulting coordinates and velocities were used for the production run of 6 ps duration with the same reaction coordinate and force constant. The reaction coordinate was then increased, the system with coordinates and velocities from the previous run was equilibrated for 2 ps followed by a 6 ps production run. This procedure was followed for the rest of the windows.

The solution simulations (both in the presence and absence of Mg²⁺) were done in an 11 Å radius sphere of 179 water molecules with stochastic boundary conditions at 300 K. The TIP3P model [33] was employed for the waters and &Aringqvist parameters were used for Mg²⁺ [34].

The number of windows used for the solution simulation in the absence of the divalent cation is 20. The starting structure is the same as in the gas phase. The water was first energy minimised for 200 steps with the reaction system fixed. Then the whole system was equilibrated at 300 K for 5 ps with the reaction coordinate set to 1 followed by a production run of 5 ps. The coordinates and velocities from previous runs were used for successive equilibration and production runs at successive windows with increasing reaction coordinate.

The simulation with Mg²⁺ ion was essentially similar. The Mg²⁺ ion was placed at the bisector of the O_ax1-P-O angle and at a distance of around 4 Å from the phosphorus atom before immersing the reaction system in the solvent sphere.

The time spans used for equilibration and production runs in each window was seen to be sufficient for our small system as revealed by checking the constancy of temperature and various energy components with time. These time scales may also be compared with those employed by other authors for similar purpose. Bash et al [35] employed 0.5 ps equilibration and production runs for a hybrid QM/MM free energy perturbation calculation involving Cl^- + CH₃Cl. Ho et al [36] showed that a 10 ps production run gave essentially the same result as a 5 ps run in their investigation of proton and hydride transfer reactions in solution.

Fig.3 Snapshot of a simulation step

Results & Discussion

The results of the umbrella sampling procedure are given in figure 4.

Fig.4 Reaction free energy profiles (in kcal/mol) as a function of reaction coordinate.

The energy profiles reveal some notable features:

There is a deep minimum around the value of 1 for the reaction coordinate implying the existence of a stable intermediate. Solvation with water molecules increases the well depth and the stability by 7.9 kcal/mol. On the otherhand, Mg²⁺ almost eliminates the intermediate and also shifts the position of the minimum towards the side on which it resides.

The barrier heights from the dissociated state are also strongly dependent on the environment. Solvation alone helps to reduce the barrier by 15 kcal/mol whereas in the presence of Mg²⁺, it goes down by 27 kcal/mol from the gas phase value and 12 kcal/mol from that of the pure solvent. On the whole, the reaction free energy profiles depicted above show very clearly the facilitation of the bond breakage process in the presence of Mg²⁺.

Reduction in the energy barrier during solvation is understandable if one considers that both of the products resulting from the disintegration of the intermediate are negatively charged. The water around the reactive complex effectively screens the repulsive coulomb interaction, resulting in a reduced energy cost. The same water environment is also responsible for the stabilization of the intermediate relative to the gas phase. A more problematic aspect is the role of Mg²⁺. Although this study clearly shows that its eletrostatic effect is sufficient to cause a major reduction in the barrier height (confirming the hypothesis put forward by Warshel [37], that electrostatic interactions are the dominant factors in enzymic catalysis), it is somewhat unclear why there should be a destabilization of the intermediate in the presence of Mg²⁺.

It is also of interest to look at the positional fluctuation of the Mg²⁺ ion as the simulation progresses. Since in the simulation we did not fix its position or form any covalent linkages, it is free to move about restrained only by Coulomb and steric interactions. In figure 5 we have plotted the drift of the ion as its distance from the central phosphorus. As the distance between the leaving group oxygen and phosphorus increases, the Mg²⁺ is also seen to be a little bit dragged along with the oxygen.

Fig.5 Distances of the leaving group oxygen and Mg²⁺ ion from the central phosphorus (in Å) plotted against time (in ps).

It is instructive to compare the results of our calculations with previously performed ab initio results for similar oxyphosphorane compounds. The problem of the existence of the intermediate is a multifaceted one. On one hand, it shows a basis set dependence. Using minimal basis set reveals the existence of a substantial minimum which disappears on using a higher level theory for the same system [7]. On the otherhand, the stability of the intermediate is also dependent on the effectiveness of charge delocalization. Replacing hydroxyl groups by methoxyl groups increases the stability. Stability is also seen to increase if instead of employing a dianionic species, one takes a monoanionic or neutral one. Solvating the reactive complex with a few water molecules is also seen to induce a stable intermediate. which has been again ascribed to charge delocalization by the water molecules [9]. Interestingly, our results show that water can stabilize the intermediate without reducing the charge on the complex, presumably through dipolar reorientation.

To be able to extend the calculation to the actual enzyme system, one needs to investigate a few more features. First of all, one needs to know the accuracy of the semi-empirical hamiltonian in modelling the reaction. Secondly, one also has to take into consideration the fact, that the force field parameters used in the molecular mechanics part may not reproduce the interaction between the classical and quantum part very well. Recently, a method for accurate calibration of the semi-empirical hamiltonian with respect to rigorous quantum calculations for a particular type of reaction was proposed and applied to study proton and hydride transfer reactions using the hybrid QM/MM free energy calculation [36,38]. We are at present working with a similar method for selecting and calibrating the most suitable semi-empirical hamiltonian.

Conclusion

In conclusion, two features of the P-O_ax bond breakage reaction of a model oxyphosphorane compound have been identified by running hybrid quantum and classical dynamics. It is found that the gas phase reaction barrier is lowered in the presence of solvent and is drastically reduced in the presence of Mg²⁺. It is also observed that the stability of the intermediate increases in pure solvent but in the presence of Mg²⁺ on the contrary it is reduced. All these observations are in qualitative agreement with previous ab initio results. Interestingly, it appears that treating the solvent and the cation classically is sufficient to produce the observed effects.

Acknowledgements

We would like to thank Dr. Leif Eriksson and Dr. Srikanta Sen for useful discussions and the Wenner-Gren Foundation for financial support to one of the authors (A.L.).

References