Computational Modeling of the Catalytic Reaction in Triosephosphate Isomerase

Abstract

We present a comprehensive analysis of the catalytic cycle of the enzyme triosephosphate isomerase (TIM), including both the reactive chemistry and the catalytic loop and side-chain motions. Combining accurate mixed quantum mechanics/molecular mechanics (QM/MM) and protein structure prediction methods, we have modeled both the structural and chemical aspects of the reversible isomerization of dihydroxyacetone phosphate (DHAP) to d-glyceraldehyde 3-phosphate (GAP), for which there is a wealth of experimental data. The conjunction of this novel computational approach with the use of the recent near-atomic resolution TIM-DHAP Michaelis complex PDB structure, 1NEY.pdb, has enabled us to obtain robust qualitative and, where available, quantitative agreement with a wide range of experimental data. Among the principal conclusions that we are able to draw are the importance of the monoanionic (as opposed to dianioic) form of the substrate phosphate group in the catalytic cycle, detailed positioning and energetics of the key catalytic residues in the active-site, the flexible nature of Glu165, which favors its direct involvement in the formation of the enediol intermediate, energetics of the open and closed form of the catalytic loop region in the presence and absence of substrate, and quantitative reproduction of various experimentally measured reaction rates, typically to within ∼1 kcal/mol. Our results are consistent with the available experimental data, and provide an initial picture as to why loop opening when GAP is the product has a higher barrier than when DHAP is the product.

Introduction

Triosephosphate isomerase (TIM) is one of the most extensively characterized enzymes in the chemical literature. TIM catalyses the reversible isomerization of dihydroxyacetone phosphate (DHAP) to d-glyceraldehyde 3-phosphate (GAP) with exceptionally high efficiency, while suppressing elimination of orthophosphate. To accomplish this isomerization, TIM first extracts a pro-R hydrogen from the C1 of DHAP and then stereospecifically introduces a proton onto the C2 atom (Figure 1). A wealth of thermodynamic, spectroscopic, kinetic, and structural studies have been carried out, providing extensive experimental data with which to compare theoretical calculations. Furthermore, the well-known large-scale catalytic loop motion in TIM (loop 6, residues 166–177), which alternates between open and closed positions, poses a substantial challenge to computational modeling efforts, one that is different from modeling enzyme reactions with relatively small-scale structural changes in the protein over the course of the catalytic cycle. The closed form of the loop is observed whenever the enzyme is ligated, and is responsible for properly aligning active-site residues fundamental to the catalytic process and for the stabilization of reaction intermediates.¹

The conversion of DHAP to GAP, a thermodynamically unfavorable but biologically significant reaction, occurs with a catalytic rate $\text{k}{}_{\mathrm{<mtext>cat</mtext>}}\text{=7.5(±0.2×10}{}^2\text{s}{}^{\mathrm{-1}}$ and a Michaelis constant $\text{K}{}_{\mathrm{<mtext>m</mtext>}}\text{=1.4(±0.1)}\text{mM}\text{.}$ The rate-limiting step appears to be a slow conformational change as described by a variety of isotope data.² In the reverse direction, with GAP as a substrate, k_cat=8.7(±0.3)×10³ s⁻¹, $\text{K}{}_{\mathrm{<mtext>m</mtext>}}\text{=0.055(±0.004)}\text{mM}\text{,}$ and the rate-limiting step appears to be a proton transfer.3., 4. Therefore, the opening of the loop on DHAP appears to be faster than 8×10³ s⁻¹, whereas on GAP it must be comparable to 7×10² s⁻¹. Recently, these kinetic measurements have been complemented with a variety of spectroscopic kinetic probes of the loop-opening rates in the presence of a variety of ligands.5., 6., 7. Information on the chemical mechanism and the key players at the active-site has been provided from a variety of enzymological, mutagenic and spectroscopic studies. Mutational changes of Lys12⁸ and Glu165⁹ have important effects in both K_m and k_cat, whereas the His95Q mutant reduces the catalytic rate by 200-fold.10., 11. A large collection of crystal structures, including the recent TIM-DHAP Michaelis complex at near-atomic resolution, 1NEY.pdb,¹² indicate a highly compact active-site. Spectroscopic and X-ray crystallographic data indicate that the substrate is primarily in the form of the ketone DHAP in situ. Short-range interactions to all three key residues, Glu165, His95 and Lys12, are seen in the ground state or Michaelis complex. It is interesting, based on this high-resolution picture of the ground state, to ask how stabilization is provided for the key intermediates. This issue is best addressed through computational studies, and must combine predictions of quantum chemical pathways with predictions of conformational freedom for key residues.

Several theoretical studies, using semiempirical and DFT-based QM/MM methods to study TIM catalysis, have provided an initial atomic-level picture of the various reaction steps.13., 14., 15., 16., 17. However, these studies were not able to utilize the recent TIM-DHAP crystal structure discussed above as a starting point. The resulting structural uncertainties, in conjunction with issues concerning the accuracy of the computational technology employed in the calculations, made it difficult to confidently draw quantitative conclusions. Furthermore, no investigation of the catalytic side-chains or loop motion was carried out. Thus, a complete, robust theoretical picture of TIM catalysis has not been produced.

Here, we use a combination of QM/MM and protein structure modeling techniques to provide a comprehensive, and quantitatively reliable, description of the TIM catalytic cycle. Our QM/MM methods enable routine treatment of quantum regions of the active-site as large as 150–200 atoms, for which full DFT geometry optimizations can be carried out with sufficiently large basis sets to guarantee accurate structures. The QM/MM interface parameterization that we employ has been tested extensively in a wide range of prior calculations on small model systems and proteins;18., 19. we have demonstrated that the errors arising from the interface are less than 1 kcal/mol $\text{(1}\text{cal}\text{=4.184}\text{J}\text{).}$ Furthermore, the overall accuracy in calculating free energies of intermediates and transition states for enzymatic reactions appears to be of the order of 2–3 kcal/mol. This technology represents a significant improvement in accuracy when compared with models and methods used in previous theoretical studies, thus enabling reliable quantitative conclusions to be drawn concerning barrier heights, protonation states, and other details of the reaction. As discussed below, the protonation state of the phosphate group on the ligand is of particular importance, and we have invested considerable effort in considering the alternatives and assessing them according to total energies obtained.

We investigate the catalytic loop motion in TIM using methods for loop and side-chain prediction,20., 21., 22., 23. based on a protein molecular mechanics force-field and a surface Generalized Born (SGB) continuum solvation model that we have developed several years ago.²⁴ The methods have been validated through an extensive series of tests in predicting loop and side-chain structures in the PDB, and have been shown to yield a remarkably high level of accuracy, even for relatively large loops (the catalytic TIM loop is 11 residues, which is within the validation range of our test suite). QM/MM calculations are used to generate charges for the ligand, which is likely to be a substantial improvement over using generic force-field parameters, particularly for the charged phosphate group.

The structural prediction methods are able to successfully generate open and closed forms of the loop, with and without the ligand present, with thermodynamics that are consistent with the entire range of experimental data. These calculations provide an initial atomic level picture of the loop opening with the ligand present, which has been proposed as the rate-determining step in the TIM catalytic cycle.

The specific reaction catalyzed by TIM that we have chosen to study is given below in Equation (1):

$$\text{E}\text{+}\text{DHAP}\text{↔}\text{EDHAP}\text{↔}\text{EGAP}\text{↔}\text{E}\text{+}\text{GAP}$$

The thermodynamic stability of the product and reactant in aqueous solution has been measured experimentally, and yields a ΔG of 1.9 kcal/mol.²⁵ Thus, the products are somewhat less stable than the reactants. In vivo, the product is essential for optimal throughput in the glycolytic pathway. Depletion of product continues to drive transformation of the reactants. Without catalysis, the reaction rate in solution is ∼10⁻⁷ s⁻¹,²⁶ which would be inadequate to serve the relevant biological function. TIM accelerates the rate by about nine orders of magnitude.

The key residues involved in the catalysis are believed to be Glu165, His95, and Lys12 (Figure 1). It is widely accepted that the first step involves the transfer of a proton from C1 of DHAP to the catalytic base, Glu165 (atom notation can be followed with Figure 1, Figure 2). From this point, three different paths have been proposed and discussed by several experimental and theoretical groups. Paths A and C involve the formation of an enediol intermediate by transfer of a proton from His95 (path A) or Glu165 (path C) to the substrate O2 atom. The enediol intermediate then forms enediolate by the transfer of a proton from O1 to His95 (path A) or Glu165 (path C). Path A, the “classical” mechanism, was proposed by Knowles et al.⁴ and has been accepted widely. Path C has been supported recently by NMR experiments,²⁷ and by ab initio active-site theoretical models.¹⁷ Path B¹³ forms an enediolate intermediate directly by internal proton transfer from O1 to O2. The final catalytic step, common to all paths, involves a proton transfer from Glu165 to C2 to form the product GAP. The proposed steps for the reaction should account for central experimental evidence: the substantial isotope washout when ²H and ³H are used as substrate.²⁸

In previous theoretical studies, the enzymatic cycle was studied using the dianionic (D state) form of the ligand phosphate group. This is the thermodynamically most stable form of the ligand in neutral solution, as the second pK_a of the phosphate group is ∼6.5. However, the free energy difference between the singly (monoanionic, M state) and doubly charged species would be less than 1 kcal/mol; furthermore, partial burial of the phosphate group in the enzyme active-site would tend to favor the protonated species, as the outer shell Born energy of the dianion (substantially larger than that of a single charged species) would be diminished. Thus, a pK_a shift from the solution value is expected. In what follows, we study the above reactions using both the singly (M) and doubly (D) charged forms of the phosphate group. We did briefly investigate the neutral species as well, which produce a substantial reduction in activation barrier; however, the pK_a of this last proton is extremely low, ∼2–3, making it highly unlikely that this species contributes substantially to the reaction mechanism. Note that even if the effect of the protein on pK_a is negligibly small and the dianion substrate is the major bound complex,²⁹ the singly charged form can still dominate the optimum reaction pathway; all that is required is that the overall activation free energy of this pathway is lower than that achieved by going through the doubly charged species exclusively. In formulating this argument, we assume that protonation is diffusion-controlled, so that the kinetic barrier to forming the monoanion from the dianion at pH 7 would be significantly lower than that observed for the overall reaction.