J. Phys. Chem. B 2007, 111, 7665-7674
7665
Substrate Orientation in 4-Oxalocrotonate Tautomerase and Its Effect on QM/MM Energy Profiles Tell Tuttle* and Walter Thiel* Max-Planck-Institut fu¨r Kohlenforschung, D-45470 Mu¨lheim an der Ruhr, Germany ReceiVed: December 14, 2006; In Final Form: March 30, 2007
The tautomerization of 2-oxo-4E-hexendioate by 4-oxalocrotonate tautomerase has been studied by quantum mechanical/molecular mechanical (QM/MM) methods using three models, A-C, with different substrate orientations. The computed QM/MM energy profiles are rather different. Various energy partitioning analyses indicate the origin of these differences and the role of the active site residues for different substrate orientations. The proposed new model C is preferred over the previously used models A and B because it combines favorable substrate binding geometries with reasonable barriers and is consistent with the experimental evidence from mutation studies concerning the catalytic ability of specific residues in the binding site, especially R11′.
Introduction 4-Oxalocrotonate tautomerase (4-OT) is an enzyme that is produced by Pseudomonas putida mt-2 (a soil bacterium), which is active in a metabolic degradation pathway to convert aromatic hydrocarbons to energy.1,2 The substrate specific to 4-OT is 2-oxo-4E-hexendioate (S, Scheme 1). Once bound, the enzyme catalyzes the isomerization reaction resulting in 2-oxo-3Ehexendioate (P). Experimentally, the isomerization reaction, when catalyzed by 4-OT, is extremely facile and proceeds with a barrier of ca. 13 kcal/mol3,4 (estimated from reported kcat values using transition state theory). The enzyme has a hexameric structure (formed as a trimer of dimers), where each monomer unit is 62 amino acids in length. The relatively small size of each subunit makes it amenable to total chemical synthesis. As such, it has been the focus of several chemical mutation studies, which provide insight into the role of individual residues that affect the catalytic activity.4-7 Such understanding is helpful in the design of mutants and synthetic analogues with new catalytic abilities (e.g., through replacement of arginine residues by citrulline).8 Previous experimental work on 4-OT3,9-11 has identified the terminal proline residue (Pro-1) as a general base (pKa ∼ 6.4 in the hydrophobic active site) which transfers the proton from C3 to C5 (for atom numbering see Scheme 1). An X-ray study12 has located three arginine residues (Arg11′, Arg39′′, and Arg61′) in the binding site that are expected to interact strongly with the dianionic substrate, and a fourth residue (Leu8′) has recently also been suggested to be catalytically important13 (primed and double-primed notation indicates different monomer chains of the enzyme). The possibility of one of the arginine residues acting as a general acid in an “acid-base catalysis” scenario has been excluded based on the experimentally determined pKa values of the three residues.14 The role of individual residues during binding and catalysis can be investigated computationally. Quantum mechanical/ molecular mechanical (QM/MM) methods are ideally suited for this purpose, and several QM/MM studies on 4-OT have already been published which address this issue.7,13,15,16 The outcome * Corresponding authors. E-mails:
[email protected] (W.T.) and
[email protected] (T.T.).
SCHEME 1: Atom Numbering of the Substrate (S) and the Isomerization Product (P)
of such QM/MM studies is influenced by many choices that need to be made with regard to computational methodology and system setup. In this paper, we examine the effect of these choices on the computed QM/MM energy profiles for 4-OT. We use three models with different substrate orientation and investigate the role of the binding-site residues in these models. Computational Methods Model Setup. The two structural models used in the initial stages of this work have been described previously.15,16 Model A was constructed for the purpose of studying the effect of mutations in the binding site on the catalytic proficiency of 4-OT,16 whereas model B was obtained from Cisneros and Yang.15 There are two main differences in the setup of models A and B: (1) the orientation of the substrate in the active site and (2) the protonation states of the histidine residues in the enzyme. In both models, the structure of the enzyme was taken from the protein data bank (PDB ID: 1BJP). This structure of 4-OT is an inactive form that contains a covalently linked inhibitor 2-oxo-3-pentynoate.12,17 The six substrates (S) were placed in the six active sites of the enzyme using the position of the covalently linked inhibitors for an initial alignment. The 6-carboxylate group (see Scheme 1 for numbering) is not present in the inhibitor, and as such its orientation is ambiguous. In model A, the 6-carboxylate group was positioned such that it could interact with Arg-11′ through H-bonding (Scheme 2a), whereas in model B a linear form of the substrate was chosen, which lacks the interaction with Arg-11′ (Scheme 2b). In model A, the protonation states of ionizable residues were assigned based on a previous theoretical study18 and crosschecked with the empirical pKa prediction program, propKa.19
10.1021/jp0685986 CCC: $37.00 © 2007 American Chemical Society Published on Web 06/14/2007
7666 J. Phys. Chem. B, Vol. 111, No. 26, 2007
Tuttle and Thiel
SCHEME 2: Alignment of the Substrate in the Binding Site: (a) Model A; (b) Model B
On the basis of this evaluation, the histidine residues were assigned a neutral total charge. There are 12 such residues in total (two per monomer) where those closest to the active site belong to a separate monomer (i.e., H6′ and H49′ are ca. 11 and 14 Å away from the active site, respectively).20 The chosen assignment leads to a total charge of -1e for each monomer (i.e., the total system charge is -6e). The addition of the six dianionic substrates results in a total system charge of -18e. In model B, all histidine residues were protonated, yielding a total system charge of -6e. In both setup procedures one of the six possible binding sites was hydrated using the droplet model centered at the C3 atom of S1. Following the hydration, the systems were equilibrated in the reactant state of the complex at the MM level of theory. For full details of the parameter development, hydration procedures, geometry optimization, and MD simulations, the reader is referred to the previous publications.15,16 In the course of this work, a third model C was constructed by a QM/MM re-equilibration of model A starting from the optimized structure of the reaction intermediate. This led to a reorientation of the Arg39′′ residue through a 180° rotation around the Cδ-N bond accompanied by a reorganization of the H-bond network in this region (see below). Models A and C are analogous in all other aspects. QM/MM Setup. The original QM/MM study by Cisneros and Yang on the model B system has been described elsewhere.15 In the current work we have used our own QM/MM methods, and in order to eliminate differences between the two approaches all model B structures were recalculated using their coordinates. The QM/MM calculations were carried out at the (BLYP21,22/TZVPP23)/CHARMM24 level of theory, unless otherwise stated. The modular program package ChemShell25,26 was employed for all QM/MM calculations, where the QM energy and gradients were provided by Turbomole27-31 (QM ) BLYP). ChemShell’s internal force field driver using the CHARMM parameter and topology data provided the MM energy and gradients. The QM/MM coupling was performed at an electrostatic embedding level32 and the boundary region was treated via a link atom approach with the charge-shift scheme.26,33 The QM region contains the substrate (S1) and the terminal proline (P1) adjacent to S1, resulting in 30 QM atoms (31 atoms including the link atom). The cut between the QM and MM regions coincides with the CHARMM charge-group definitions for the P1 residue such that the QM and MM regions both have integer charge. For the QM/MM geometry optimizations a subset of the total system is optimized while the remainder of the atomic positions are frozen at the equilibrated geometry. The active subset for all optimizations was defined from the initial reactant complex geometry using a distance criterion, whereby any residue that contains an atom within 10.0 Å of any atom of S1 is included in the active subset. The resulting active region contains ca. 1100 atoms (1114 for model A), about one tenth of the total system size (ca. 13 500 atoms).
Results and Discussion Structural Comparison of Model A and Model B. The primary structural difference between model A and model B is the orientation of the substrate in the binding site (Figures 1 and 2, respectively). In model A (Figure 1) there are two strong H-bonds formed between the 6-carboxylate of the substrate and R11′. In addition, the 6-carboxylate terminal also interacts with the L8′ residue of the same chain. However, in model B (Figure 2) the interaction with R11′ is missing; the closest contact between the residue and the substrate is ca. 3.9 Å. The interaction with R39′′ is also strongly altered in the different orientations. In model A R39′′ interacts directly with the 2-carbonyl O atom of the substrate, and in the reactant and product structures an H-bond is formed to one of the 1-carboxylate O atoms. In model B there is no interaction between R39′′ and the 2-carbonyl O atom. Rather, the 2-oxo atom in this system is weakly bound to a water molecule (see Scheme 2). Thus, in model B, R39′′ binds exclusively to the 1-carboxylate O atoms of the substrate. In both setups there is a clear H-bonding interaction between the 1-carboxylate O atoms and R61′. In model A this interaction involves both O atoms, forming two strong H-bonds. In model B only one of the O atoms is bound to R61′ (the other is bound to R39′′) although the interaction is still strong. The stronger pairing between the terminal carboxylate groups of the substrate and the arginine residues in the binding site suggests that the substrate will be more strongly bound in model A than in model B. Methodological Comparisons. The relative energies published previously for model B differ significantly from those obtained with model A. The reaction profile for the isomerization reaction is shown in Figure 3, where the potential energy values (∆E) obtained in the current work for model A and in the previous work on model B13 are given. In the model A orientation with strongly bound substrate, the calculated barriers are significantly higher than the activation energy estimated from experiment (ca. 13 kcal/mol). On the other hand, the model B system with its less stable binding site has a barrier that is comparable to experiment. To analyze this difference further we consider the different factors that contribute to the calculated reaction energy profiles. The QM/MM setup used in this work differs from that used in the previous model B study13 in several ways: (1) Both the QM and MM levels of theory are different - this work: QM ) RI-BLYP/TZVPP; MM ) CHARMM (with unique parametrization for terminal proline and S); previous work: QM ) HF/3-21G (with single point calculations at B3LYP/6-31G*); MM ) AMBER (with unique parametrization for S). (2) The coupling schemes are different - this work: H-link atom approach with charge scaling; previous work: pseudobond model. (3) The optimization algorithms are different.
Substrate Orientation in 4-OT
J. Phys. Chem. B, Vol. 111, No. 26, 2007 7667
Figure 1. Model A optimized binding site structures. Distances in Å. The substrate is in ball-and-stick representation while the surrounding residues are in Licorice representation. Color code for atoms: H, white; C, cyan; N, blue; O, red. Noncovalent interactions are indicated by black lines.
The effects of the first two factors were tested by performing single point calculations on the model B structures at our QM/ MM level of theory. These single-point calculations yield energies that compare well (within 1.5 kcal/mol) to the B3LYP34 values reported (Table 1). This implies that the differences in the QM/MM methods (i.e., level of theory and coupling schemes) have only a minor effect on the results. As a test of the optimization algorithms used in locating the transition states we re-optimized the first stage of the isomerization reaction (i.e., react f INT) for the model B structures. This process results in a slight (1 kcal/mol) lowering of the barrier (Table 1, BLYP-opt column) but did not cause any major structural rearrangements in the active site. Therefore, for a given system, both algorithms (HDLC_Opt in ChemShell and the nudged elastic band method used by Yang et al.) converge to
similar results. Thus, we can exclude the optimization procedure as a major source of error. Decomposition of QM/MM Energies. The QM/MM framework that we employ involves an electronic embedding scheme. The QM/MM energy is composed of three contributions:
ETotal ) EQM + EMM + EQM/MM
(1)
where ETotal is the total energy of the full QM/MM system, EQM is the energy of the QM subsystem; EMM is the energy of the MM subsystem, and EQM/MM contains the interactions between QM and MM subsystems (electrostatic, vdW and bonded interactions). However, in practice, the electrostatic interaction between the QM and MM subsystems is included in the QM calculation, through the addition of the MM point charges into the QM Hamiltonian, whereas the bonding and van der Waals
7668 J. Phys. Chem. B, Vol. 111, No. 26, 2007
Tuttle and Thiel
Figure 2. Model B optimized binding site structures. Distances in Å. The substrate is in ball-and-stick representation while the surrounding residues are in Licorice representation. Color code for atoms: H, white; C, cyan; N, blue; O, red. Noncovalent interactions are indicated by black lines.
(vdW) interactions between the QM and MM subsystems are determined in the MM calculation. Thus, for practical purposes the total energy may be defined as
ETotal ) E(QM,MM) + E(MM,QM)
(2)
where E(QM,MM) is the sum of EQM and the energy resulting from the electrostatic interaction between the QM and MM subsystems and E(MM,QM) is the sum of EMM and the vdW and bonded interactions between the MM and QM subsystems. A decomposition analysis according to eqs 1 and 2 shows that in both models there is a substantial contribution from the destabilization of the MM system (∆EMM) to the energy of the TS and intermediate structures. Indeed, in the model B structures, the QM energies that include the electrostatic interaction with the environment (∆E(QM,MM)), would yield a qualitatively different reaction profile, in which the intermediate
is stabilized below the reactant and the barrier to the first step in the reaction (TS1) is less than 1 kcal/mol. Clearly, the QM/ MM interactions (∆EQM/MM) give a major contribution in both models which mostly comes from the electrostatic interactions (i.e., ∆E(QM,MM) - ∆EQM). The contribution of the other QM/ MM interactions (∆E(MM,QM) - ∆EMM) is more pronounced in model B than in the model A system, mostly due to differences in the vdW QM/MM interactions (see Supporting Information for details). Effect of Individual Residues in the Binding Site on the Relative Energies. The strong electrostatic interaction between the MM region and the QM region was investigated on a perresidue basis by charge deletion analysis, i.e., by comparing between E(QM,MM) of the full system and E(QM,MM) of the system with the charges on the selected residue deleted. Their difference (∆E(QM,MM)) is a measure for the interaction between the residue
Substrate Orientation in 4-OT
J. Phys. Chem. B, Vol. 111, No. 26, 2007 7669 TABLE 3: The Effect of Individual Residues (R ) Arginine, L ) Leucine; See Figures 1 and 2)a react
TS1
INT
TS2
Prod
residue ∆E(QM,MM) ∆∆E(QM,MM) ∆∆E(QM,MM) ∆∆E(QM,MM) ∆∆E(QM,MM)
Figure 3. Reaction profile for the isomerization reaction of S in 4-OT. Values in bold font are for model A. Other values correspond to model B energies from a previous study.13 All values in kcal/mol.
TABLE 1: Relative Energies of the Model B Structures at Different QM/MM Levels of Theorya react TS1 INT TS2 Prod
HF
B3LYP
BLYP
BLYP-opt
0 16.9 9.3 20.6 1.0
0 14.3 11.6 17.7 1.8
0 12.9 10.7 18.4 2.0
0 11.9 6.7
a Energies in kcal/mol, relative to the reactant state. HF and B3LYP values were reported previously.13
TABLE 2: Decomposition of ∆ETotal for the Calculated Reaction Patha structure ∆ETotal ∆E(QM,MM) ∆E(MM,QM) ∆EQM ∆EMM ∆EQM/MM react TS1 INT TS2 prod
0.0 29.1 26.5 29.5 -8.8
0.0 18.3 12.3 18.5 -8.6
Model A 0.0 10.8 14.1 11.1 -0.2
0.0 52.3 58.0 46.8 -3.1
0.0 8.8 9.9 7.0 -1.4
0.0 -32.1 -41.4 -24.2 -4.3
react TS1 INT TS2 prod
0.0 12.9 10.7 18.4 2.0
0.0 0.3 -10.8 3.0 -2.1
Model B 0.0 12.7 21.5 15.5 4.1
0.0 32.0 45.7 36.0 1.7
0.0 6.4 12.9 12.0 3.5
0.0 -25.5 -47.9 -26.0 -3.2
a
Energies in kcal/mol, relative to the reactant state.
and the QM region (i.e., both the direct interaction energy and the work done in polarizing the QM density). Table 3 lists the ∆E(QM,MM) values obtained for the reactant and their variation during the reaction. For all arginine residues, the ∆E(QM,MM) values are large and negative, as can be expected from the interaction between the positive arginine residues and the negative QM region, while those for the neutral leucine residue (L8′) are about an order of magnitude smaller. The ∆E(QM,MM) results for model A and model B are qualitatively similar, except for the contribution of R11′, which is twice as large in model A. This is a direct consequence of the different orientations of the substrate in the binding site (see Figures 1 and 2). Concerning the variation of the interaction energies during the reaction, the most striking difference between the two models is the destabilizing effect of R11′ in model A for all stationary points, relative to the reactant. In model B, R11′ has only a minimal (
