Hybrid QM/MM Studies of the Reaction Mechanism of Dihydrofolate Reductase

While much has been learned about enzyme mechanisms at the molecular level from x-ray crystallography, the picture is largely static, while other experimental techniques offer a more thermodynamic view but only at the macroscopic level. In principle, computer simulation provides a means to complement experimental data and provide an accurate energetic representation of the chemical events at the reaction centre, and of the energetic and dynamical roles of the enzyme and its components. In particular, it offers a means to study the "invisible" players, protons and water, not seen by crystallography. All these players and properties may play a role in enzyme catalysis: the problem is how to decipher for a particular enzyme of interest. Delineation of the key contributions for a molecular system of typical enzyme size (i.e. at least 3,000 protein atoms plus, at least, an equivalent number of solvent molecules) is beyond the feasibility of any current computational method; no current method can provide this reliably at the level now routinely achievable in small-molecule chemistry. Rather, current methods offer "windows" which can be divided into three main classes (QM level and type; MM/MD; and hybrid QM/MM+MD) depending on the quality of the electronic description of the molecule; the quality of the representation of its chemical environment; and the degree of dynamic sampling. We are pioneering a multi-method approach using all these methods to study the reaction mechanism of the enzyme dihydrofolate reductase (DHFR), and, by comparative evaluation, to define the ranges of applicability of the different methods to large complex protein systems. In addition to scientific interest and applied interest as a major drug target, DHFR is well recognized as a test system. This is because the catalytic mechanism and ligand-binding properties have proven so refractory to experimental and computational investigation over 30-40 years, that there is an expectation their elucidation will offer general insights into enzyme mechanisms and how to simulate them.

Principal Investigator

Dr Jill E. Gready
Computational Proteomics Group
John Curtin School of Medical Research



Facilities Used



Dr Peter L. Cummins
Dr Stephen P. Greatbanks
Computational Proteomics Group
John Curtin School of Medical Research

Dr Ivan Rostov
ANU Supercomputer Facility

Dr Alistair P. Rendell
Dept of Computer Science

Prof Kritsana Sagarik
School of Chemistry
Suranaree University of Technology

RFCD Codes

250106, 250601, 250699, 270108

Significant Achievements, Anticipated Outcomes and Future Work

Background. DHFR catalyses the reduction of folate to dihydrofolate(DHF) and tetrahydrofolate (THF) via a two-step reaction of protonation, and hydride-ion transfer from the cofactor NADPH. Key aspects of the mechanism have long remained elusive, specifically: whether the carboxylate group of the conserved active-site Asp or Glu residue is protonated or ionized during the reaction; the source of the protons; and the role of active-site water molecules (especially W206, E. coli DHFR numbering). These issues have been difficult to define computationally because of DHFR complexities: it seems to involve many, if not all, protein complexities (proton shifts, active-site water, loop-closing mechanisms and dynamics, highly non-isotropic global and active-site electrostatics); it is an old essential and ubiquitous enzyme with an odd chemically-taxing role, catalysing two hydride-ion transfer reactions on two different carbons of a complex N-heterobicyclic p-system (pterin ring) in two oxidation states, and with an attendant requirement for protonations of two different ring nitrogens; and, finally, accurate electronic description of these large ring systems and the active moiety (nicotinamide ring) of the cofactor NADPH, and changes during the reaction, is non-trivial. There is no basis for assuming these effects are additive, so the challenge is to define how they interact, and to calculate their co-operative quantitative contributions.

Results. With respect to the main questions and complexities, we have established so far that:

(i) Polarization by DHFR of its substrates is not sufficient to favour hydride-ion transfer at C7 (folate) or C6 (DHF) without pre-protonation, as reported in the literature from earlier DFT calculations. The studies used both point-charge embedded ab initio QM (SCF, MP2 and various DFT models) and linear-scaling semiempirical QM (MOZYME) calculations [1,2]. Large substrate polarizations from the DFT results were shown to be artifacts of the DFT description of large anions, i.e. they represent a clear example of method breakdown. However, a new finding was that polarization at the active carbon (C4) of the nicotinamide ring does activate C4 for hydride-ion transfer. The MOZYME calculations [2], also allowed, uniquely, evaluation of charge-transfer effects to the enzyme; these were small in agreement with other MOZYME studies on the reaction path (see report for x11).

(ii) A combined ab initio QM and QM/MM+MD study [3] of the energetics of the H-bond interaction between the active-site carboxylate (i.e. ionized form) and a protonated substrate (a folate mimic, 8-methyl-pterin) predicted that the neutral-pair form (i.e. proton transferred from N3 of substrate to form the carboxylic acid) is more stable than the ion-pair form, the form universally assumed in DHFR discussions in the literature. This is the case both in vacuum and in the enzyme although this relative stability is reduced to ~5 kcal/mol by the more favourable electrostatic interactions made by the ion pair with the enzyme. Such a neutral-pair form does not activate the pterin ring for hydride-ion transfer, as required by the pre-protonation hypothesis for catalysis. However, the results were suggestive of a small free energy gap and the possibility of a low energy barrier for the proton transfer, a mechanism proposed in the literature as operative in many enzyme reactions. Followup ab initio SCRF calculations [4] indeed established a low-barrier hydrogen bond (LBHB): the ion pair corresponded to a shallow local minimum at the transition state (TS) between the two neutral-pair complexes. Thus, the H-bonded system has maximum probability of being found in the ion-pair form. As analysis of charge distributions indicates the TS is highly cationic, a mechanism for activation of substrate towards hydride-ion transfer similar to that of the formal cation seemed possible.

(iii) However, using a combination of methods as for (ii), extension of these studies to consider both the role of protonation and a conserved active-site water molecule for the binding and reaction of folate and dihydrofolate provided another picture [6,7]. Results showed that a conserved water molecule (W206, E. coli DHFR) that is H bonded to both the OD2 oxygen atom of the carboxyl (Asp) sidechain and O4 of the pterin/dihydropterin ring, is critically important, and appears to determine preferred protonation sites for the enzyme-bound substrates [6]. Single protonation of the system does not lead to the N8 (folate) and N5 (DHF) protonated forms required for a pre-protonation activation mechanism, but these forms are the most stable after addition of two protons; the other proton resides on the carboxyl group. This result also appears to have indirect experimental support: calculated H-bond distances for W206 and protonated OD2 complexes (OD2–N3 distances < OD2–NA2 distances) correlate well with those observed in many x-ray structures. Analogous studies for the reaction at the enthalpic level (ab initio QM) showed that reduction takes place only if the active-site Asp27 is protonated, but that at the free-energy level (QM/MM +MD) the catalytic advantage of a protonated Asp27 residue is small [7]. For both ionized and neutral states of Asp27, analysis of the electrostatic interaction energies between the QM and MM regions in the MD simulations showed that the enzymic environment (MM region) plays a far greater role in stabilizing the final products than in stabilizing the TS complex. However, although for the larger QM region model (85 atoms), the range of activation free energies (10 - 12 kcal/mol) is similar to an experimental value (13 kcal/mol), in general quantitative agreement with experimental free energies was poor.

In summary using our mixed-method approach, we have established that a balanced description of the both the active-site and wider-enzyme environment needs to be built up. Provisionally we have established that: the conserved active-site water (W206, E. coli DHFR) is critically important for correct protonation of the enzyme-bound substrates, but plays little role in the actual hydride-ion step; Asp/Glu (Asp27, E. coli DHFR) seems to be needed to be protonated for both steps; and the enzymic environment plays a greater role in stabilizing the final products than the TS complex. The last result is novel and interesting as the reaction is rate-limited by product dissociation and is functionally irreversible. However, even the most recent calculations do not allow firm conclusions to be drawn.

Future work. Work along the following lines will be undertaken to test these results, to quantitate energetic components of the two reactions, and also to further develop and evaluate methods (see below and work in conjunction with project u51/d52). A large series of multiple-trajectory QM/MM calculations (200-500 ps each) will be performed on both protonation and reaction steps with the aim of testing the above provisional conclusions and the other mechanistic issues, as well as understanding the convergence of the results with respect to system size and composition. A systematic buildup of the QM region to include complete substrate and cofactor molecules as well as additional active-site residues is planned, particularly residues perpendicular to the pteridine-ring plane. This will allow estimation of long-range electrostatic effects from the cofactor, and possible p-p stacking interactions. As this would require QM regions of up to 250 atoms with 9 link atoms, concepts being developed in project x11 for choosing QM boundaries will be employed. These studies will be complemented by very large ab initio fragment calculations on the QM regions. Our experience (and that of others in the literature) is that geometry optimizations of such systems with "untethered" ends (i.e. lacking constraints which would apply in the enzyme) are unpredictable; this will likely to much worse for larger QM-region analogues (i.e. up to 9 broken bonds). To date we have used fixed-atom constraints, but these are very artificial: we are now exploring the use of methods such as ONIOM in GAUSSIAN98, which allow imposition of softer enzyme-like potentials to restrict the conformational space of the fragment complexes.


Computational Techniques Used

Methods and programs. Three main methodologies are being used: (i) SE-QM/MM with MD simulation as implemented in our locally-developed MOPS program (see u51/d52). This includes recent parallelization under MPI in the APAC National Facility Compaq SC version of all compute-intensive terms except the QM SCF, and for running of multiple trajectories to improve the efficiency of free energy sampling. Other enhancements under development, which will be used and tested in this project, include: corrections to the SE-QM force fields to prevent both H-bond and covalent bond breaking, a significant problem for large QM regions; and use of system-specific parameters to recalibrate the SE-QM method for the specific enzymic reaction, to correct for the often large errors in reaction energies and activation energies from SE-QM. (ii) Linear-scaling MOZYME within MOPAC for all-enzyme QM calculations, as being developed and tested in project x11. (iii) GAUSSIAN98 for ab initio QM (including with implicit solvent methods, e.g. SCRF) and ONIOM (a type of QM/MM method) calculations, which was tuned and maintained for the VPP300 under the ANU-Fujitsu Area 3 contract, and is maintained for the APAC National Facility Compaq SC by National Facility staff.

Capabilities. Typical size-dependent limitations for these methods with current supercomputer resources are: 50-100 atoms for ab initio QM (e.g. SCF, DFT) calculations with full geometry optimization, including with an approximate environment (point-charge, implicit solvent); several hundred atoms for conventional semiempirical QM (e.g. AM1, PM3) calculations in MOPAC with full geometry optimization, including with implicit solvent (COSMO) model; 10,000 atoms or more for linear-scaling (MOZYME) SE-QM calculations with very limited geometry optimization, but including implicit solvent (COSMO) model; 30,000 atoms or more (enzyme + solvent) for multiple trajectory MD simulations (MOPS) of nanoseconds (ns) each; and 10,000-20,000 atoms for multiple trajectory simulations (MOPS) of ~500 ps each with a QM/MM potential (up to 200 atoms in the SE-QM region).


Publications, Awards and External Funding

1. SP Greatbanks, JE Gready, AC Limaye, AP Rendell. Enzyme polarization of substrates of dihydrofolate reductase by different theoretical methods. Proteins 37, 1999, 157-165.

2. SP Greatbanks, JE Gready, AC Limaye, AP Rendell. Comparison of enzyme polarization of ligands and charge-transfer effects for dihydrofolate reductase using point-charge embedded ab initio QM and linear-scaling semiempirical QM methods. J Comput Chem 21, 2000, 788-811.

3. PL Cummins, JE Gready. Combined quantum and molecular mechanics (QM/MM) study of the ionization state of 8-methylpterin substrate bound to dihydrofolate reductase. J Phys Chem B 104, 2000, 4503-4510.

4. PL Cummins, JE Gready. QM/MM and SCRF studies of the ionization state of 8-methylpterin substrate bound to dihydrofolate reductase: existence of a low-barrier hydrogen bond. J Mol Graph Model 18, 2000, 42-49.

5. JE Gready, PL Cummins. Mechanism-based substrates and inhibitors of dihydrofolate reductase. in Free Energy Calculations in Rational Drug Design, Reddy MR, Erion MD, eds., Kluwer Academic/Plenum Publishers, 2001, pp. 343-364.

6. PL Cummins, JE Gready. Energetically most likely substrate and active-site protonation sites and pathways in the catalytic mechanism of dihydrofolate reductase. J Am Chem Soc 123, 2001, 3418-3428.

7. PL Cummins, SP Greatbanks, AP Rendell, JE Gready. Computational methods for the study of enzymic reaction mechanisms I: application to the hydride transfer step in the catalysis of dihydrofolate reductase. J Phys Chem B, 2002, submitted.