©1996 by The American Society for Biochemistry and Molecular Biology, Inc.
Computer Simulation of the Triosephosphate Isomerase Catalyzed Reaction (*)

(Received for publication, April 10, 1995; and in revised form, January 22, 1996)

Johan Åqvist (§) Michael Fothergill (¶)

From the Department of Molecular Biology, Uppsala University, Biomedical Centre, Box 590, S-75124 Uppsala, Sweden

ABSTRACT
INTRODUCTION
MATERIALS AND METHODS
RESULTS AND DISCUSSION
FOOTNOTES
ACKNOWLEDGEMENTS
REFERENCES

ABSTRACT

A major challenge for theoretical simulation methods is the calculation of enzymic reaction rates directly from the three-dimensional protein structure together with some idea of the chemical reaction mechanism. Here, we report the evaluation of a complete free energy profile for all the elementary steps of the triosephosphate isomerase catalyzed reaction using such an approach. The results are compatible with available experimental data and also suggest which of the possible reaction intermediates is kinetically observable. In addition to previously identified catalytic residues, the simulations show that a crystallographically observed active site water molecule plays an important role during catalysis and an intersubunit interaction that could explain the low activity of the monomeric enzyme is also observed. The calculations clearly demonstrate the important catalytic effects associated with stabilization of charged high energy intermediates and reduction of reorganization energy, which are likely to be general principles of enzyme catalyzed charge transfer and separation reactions.


INTRODUCTION

A wealth of enzymological (1, 2, 3, 4) and structural (5, 6) information has been obtained over the years for the triose phosphate isomerase (TIM) (^1)catalyzed reaction, and investigations of the uncatalyzed interconversion between dihydroxyacetone phosphate (DHAP) and glyceraldehyde 3-phosphate (GAP) have also been reported(7) . The current view of the functional mechanism of the enzyme that has emerged from these studies can be summarized as follows (see Fig. 1). (i) After binding of DHAP one of the C-1 protons (1-pro-R) is abstracted by the catalytic base Glu, yielding a corresponding enediolate species. (ii) This enediolate is then protonated at O-2 by the neutral imidazole ring of His, giving a doubly protonated enediol and an imidazolate anion. (iii) His recaptures a proton from the O-1 oxygen again yielding an enediolate species, now lacking the proton at O-1. (iv) Protonation of this enediolate at C-2 by the carboxyl group of Glu, thus restoring the general base to give the product D-GAP, which then is released from the enzyme. TIM is a fast enzyme that is limited by product release at physiological substrate and product concentrations and with internal rate constants on the order of 10^4 s(3, 7) . As a comparison the rate-limiting step for nonenzymatic conversion of DHAP is about 10^5 s, in which case the mechanism involves intramolecular proton abstraction by the phosphate group(7) .


Figure 1: The catalytic mechanism of TIM is represented as a four-step process (excluding binding and release steps). Five VB configurations ((1)-(5)) are used to describe the reaction.



There are, however, still several open questions regarding the detailed mechanism that require an answer before we could claim a complete understanding of this enzyme. All of the proposed elementary steps cannot be discerned experimentally, and one instead observes only one kinetic intermediate state, but its exact nature is not known. The peculiar use of a neutral imidazole as proton donor has also attracted much attention (8, 9) and was actually a prediction from theoretical calculations of a rather unexpected result(9) . It has, however, recently been disputed in a study by Kollman and co-workers(10) . A number of other theoretical studies of TIM have been reported earlier dealing with amino acid mutations (11) and conformational changes (12) as well as quantum and molecular mechanics calculations on the catalytic reaction(13) .

Here, we report the calculation of a complete free energy profile for the proposed ``imidazolate'' reaction mechanism(4, 6, 9) using molecular dynamics (MD) free energy perturbation (FEP) simulations in combination with the empirical valence bond (EVB) method (14, 15) . Because some controversy regarding the detailed mechanism of TIM still exists, our goal here is to try to evaluate the energetics implied by the seemingly favored one in order to examine its feasibility. The x-ray coordinates of the yeast enzyme in complex with the inhibitor phosphoglycolohydroxamate (6) were used as the starting point. This inhibitor closely resembles one of the proposed (high energy) intermediates of the reaction, and no major modelling is required to construct the different reaction species (cf.Fig. 1) from it.


MATERIALS AND METHODS

The TIM reaction is described here by the EVB model(14, 15) . This model treats the reaction (or rather each reaction step) as a conversion between different valence bond (VB) states or resonance structures. Each of these pure states is modelled by a molecular mechanics force field, whereas intervening configurations are obtained by mixing of the VB states. The EVB method has the advantage that it is easily incorporated into an MD framework where the reaction path can be mapped out by the FEP procedure. Furthermore, as will be discussed below, the method allows calibration against known experimental data because the VB states have a clear physical meaning. Accordingly, thermodynamic data for model reactions in water can be used to fine tune the EVB reaction surface(15) .

In the present case, the reaction path is described by the five bonding configurations or resonance structures depicted in Fig. 1. Each of these diabatic states ((1)-(5)) is represented by an analytical force field of the standard molecular mechanics type(16, 17) . In the cases where standard force field parameters(16) , viz. partial charges, are not available, they were derived from AM1/PM3 calculations using the AMSOL program (18) and merged with the GROMOS group parameters in order to achieve consistency with them. Fig. 2. shows the charge distributions that were used for the substrate in each of the five VB states. The parameters for dipolar groups of the substrate are very close to the standard GROMOS ones. This is, of course, no guarantee that they are correct unless the corresponding solvation free energies have actually been verified(19) . However, for dipolar groups the absolute values of the solvation energies are small, and the corresponding errors will likewise be small. Moreover, there is clearly some cancellation of errors inherent in the fact that we here only consider the difference between solution and enzyme. For charged groups, on the other hand, solvation energies are large, and the nonbonded parameters are more crucial for energetics. Therefore we have employed parameters for the negatively charged oxygens that have been carefully calibrated earlier against hydration free energies(17) . The other van der Waals' parameters (16) pertaining to charged groups (e.g. for N and P) have also been checked. (^2)We are thus rather confident that the energetics associated with the present force field(s) is reasonable, and it should perhaps again be emphasized that (i) we only examine the difference between water and enzyme and (ii) the total free energy difference between different EVB resonance structures is calibrated to match experimental solution data prior to carrying out the enzyme calculations (see below). Bonds within the reacting fragments (defined as the depicted atoms in Fig. 1) were represented by Morse potentials using standard values for bond lengths and dissociation energies(20) . The CH(2)OPO(3) moiety was, however, represented with the standard harmonic bond functions (16) because these bonds are not broken during the reaction.


Figure 2: Partial atomic charges used in the calculations for the substrate in the various resonance forms. The standard GROMOS van der Waals' parameters (16) were used except for the oxygens bearing a large negative charge (q > 0.5) for which the parameters derived in (17) were used (these have been calibrated to reproduce hydration free energies for carboxylate ions).



Although the force field(s) above can describe the interactions between the fragments of each state and with the surrounding medium, as well as the energies involved in distorting geometries from their equilibrium, it does not contain information about the relative energies of the fragments in vacuum. That is to say that there is an energy difference, alpha, between any two fragments that is related to their heats of formation that is not included in the molecular mechanics force field. Furthermore, standard force fields cannot describe the quantum mechanical coupling between states reflected by off-diagonal terms of the hamiltonian. The above information can, however, be obtained either by quantum mechanical calculations or by gas phase or solution experiments. The latter alternative, to use experimental data to calibrate the above mentioned gas phase parameters of the system hamiltonian (which we refer to as EVB parameters, cf.Table 1), is the essence of the EVB method(14, 15) . In particular, it has proven useful to calibrate the potential energy surface with respect to solution experiments by simulating suitable reference reactions in water. The parameters that are subject to this calibration are the above-mentioned gas phase energy differences (alpha values) between the VB states and the off-diagonal matrix elements H (see below). The former are mainly associated with reaction free energies, whereas the latter mainly affect activation barriers. The resulting parameters are then directly transferred to simulations of the corresponding reaction steps in the solvated enzyme. The advantage with this approach is that the energetics in an aqueous environment can be reproduced exactly and one then ``only'' has to be able to model the substitution of water by protein (note, however, that the choice of reference reactions corresponds to the assumed mechanism in the enzyme and not necessarily to that of the uncatalyzed reaction, which may proceed by a different mechanism). The alternative would be to rely on a (semi-empirical) molecular orbital model to correctly reproduce the gas phase energetics and then ``substitute'' vacuum for the actual environment in the enzyme by coupling to a molecular mechanics force field(9) . This may involve more serious transferability problems, and the idea with the EVB reference reaction calibration is thus to diminish these.



In the present case some of the data needed for calibration is not directly available from experiment due to the ``high energy'' nature of the proposed intermediates, e.g. the solution pK values of hydrogens on C-1 and C-2 as well as the associated activation barriers. Starting from the pK of acetone in water of 19.2(21) , these pKs can, however, be estimated by the substituent effects of a hydroxyl group and the phosphate dianion. The latter interaction, i.e. between the enediolate and phosphate anions, has been determined to be about one pK unit in water (7) . In order to estimate the additional effect of a hydroxyl group on the pK of acetone, we carried out semi-empirical SCF calculations in a homogeneous dielectricum ( = 80) using the AMSOL program(18) . Calculations with the AM1-SM2 and PM3-SM3 hamiltonians yield a pK shift of -3.8 and -5.4 pK-units, respectively. Taking into account the observed effect of the phosphate group (7) and after correcting for the number of equivalent protons, we arrive at a pK value of 17.2 for abstracting a C1 proton from DHAP, using the AM1-SM2 result, and 15.6 using the PM3-SM3 result. The former of these values seems most reasonable in that it gives a free energy relationship that nicely fits experimental data on acetone and DHAP deprotonation from independent sources(7, 22) , as shown in Fig. 3(our estimate of the pK at 17.2 is also rather close to the value inferred by Richard(7) ). The corresponding pKvalue of 15.9 for LGAP is related by the observed equilibrium constant of to that of DHAP(7) . By interpolation in the diagram, the free energy relationship of Fig. 3can also be used to deduce the activation barriers for proton exchanges with the carboxylate moiety of glutamic acid in aqueous solution (pK = 4.1), which are used as reference reactions for the first and fourth step of the enzymic conversion (cf.Table 1). The second and third steps of the enzyme reaction, according to the mechanism of Fig. 1, involve proton transfer between the substrate oxygens and the imidazole ring of His. The pK difference between donor and acceptor in solution is about 2 pK units, with the neutral imidazole ionizing at the higher pK of 14 (9) and the enediol at 12(22) . The proton transfer barriers for these reference reactions can, in the same way as above, be determined quite accurately from a free energy relationship based on experimental data of Eigen and co-workers (Table IV of (23) and Table V of (24) ), and the corresponding energetics are listed in Table 1.


Figure 3: Free energy relationship (activation barrier versus DeltapK) comprising experimental data on acetone deprotonation by various bases ((22) , squares) and on base catalyzed phosphate elimination reactions (7) of DHAP (circles) and LGAP (triangles) using the estimated pKvalues of DHAP and LGAP. The latter pK was related to the former by the observed equilibrium constant of between LGAP and DHAP in solution(7) . All values have been corrected for the number of equivalent protons.



Hence, the energetics of the four reaction steps as they would occur in aqueous solution, assuming a nonconcerted imidazolate mechanism, can be determined as described above, and the resulting free energy diagram is shown in Fig. 4(curve a). Each of these steps is then simulated in a separate MD calculation with the reacting fragments immersed in a sphere of water. The FEP technique (25, 26, 27) is used to drive the system between the five different VB configurations (Fig. 1) while mapping out the free energy profile (potential of mean force along the reaction coordinate) on the actual ground state potential surface by the potential of mean force procedure that has been described elsewhere(14, 15, 17) . The ground state energy is obtained for any given configuration of the system (set of coordinates) by mixing the VB states and solving the corresponding secular equation (14, 15, 17) . These simulations of the uncatalyzed reference reactions are used to obtain values of the (unknown) gas phase free energy differences, alpha, between the VB states and of the off-diagonal hamiltonian matrix elements (the latter are represented by functions of the form µ = Ae


Figure 4: a, free energy diagram for the TIM reaction mechanism in aqueous solution based on experimental data (7, 21, 22, 23, 24) and semiempirical AMSOL (18) calculations. EVB/FEP/MD calibrations of each of these steps in water were carried out so that the resulting free energy profile reproduces the depicted energetics (14, 15) . b, free energy diagram for the corresponding reaction steps in the enzyme. The transition states associated with the four reaction steps are denoted (1)-(4). The top panel also shows the detailed free energy profile for the first enzymic step with the abscissa denoting the energy gap between the (1) and (2) diabatic states.



The simulation characteristics were the same as in (17) and (28) . For each step in water and in the protein, a 70-ps trajectory was calculated using about 30 values of the FEP coupling parameter , yielding a total simulation time of 560 ps excluding equilibration. At each -point the first 40% of the data was discarded. The simulations were carried out at approximately 300 K with an MD time step of 0.002 ps. Spherical systems of radius 15 Å were used with the boundaries restrained by the SCAAS model (29) (water molecules) and by harmonic restraints to crystallographic positions (protein atoms). The spheres were centered on the C-1 carbon atom both in the protein and water simulations. The substrate and charged protein groups interacted with everything within the simulation sphere, whereas cut-offs as specified in (28) were applied to other interactions. The convergence error bars obtained by forwards and backwards integration of the same trajectories were examined and found to range between 0.5 and 1.5 kcal/mol. In addition, some extra runs with different initial conditions were also carried out to try to assess the accuracy of the results, and these gave very similar error ranges. Although the above values give an estimate of the MD convergence errors (per reaction step), it should be noted that an equally important source of error for the resulting enzyme energetics is the accuracy of the solution data used for calibration. These errors are unfortunately more difficult to estimate because they are determined by the accuracy of the pK(a) values for the different species as well as the free energy relationships used to extract barrier heights ( Fig. 3and Refs. 23, 24).


RESULTS AND DISCUSSION

The resulting free energy profile from the EVB/FEP/MD calculations is summarized in the diagram of Fig. 4(curve b). The large catalytic effect of enzyme is evident, and it can be seen that all three intermediates along with their flanking transition states are stabilized relative to the reactant and product by some 15 kcal/mol compared with the uncatalyzed reactions. Of the four intervening free energy barriers, one finds that those associated with proton transfer to and from the substrate carbons limit the (internal) rate of conversion. The simulations place both (2) and (4) at lower energy than the doubly protonated species, with (2) as the most stable intermediate along the path. This would thus be the main contributor to the experimentally observed kinetic intermediate and corresponds to the case with O-1 protonated. For direct comparison with experimental results, Table 2lists the rates calculated from the transition state theory equation using a transmission factor of unity assuming that (2) is the only kinetically significant (observable) intermediate. The calculated free energy profile thus appears to be in reasonable agreement with available kinetic data(3, 7) . Both the magnitude of the kinetically significant internal activation barriers (12.5 kcal/mol) and the corresponding equilibrium constants are fairly close to the experimental estimates. The theoretical calculations thus reproduce the catalytic effect of the enzyme remarkably well, within the accuracy limits that we can estimate (see above).



The origin of the catalytic power of TIM has been the subject of much discussion(4, 6, 9, 10, 11) . We find here that there are two general effects at work: (i) stabilization of charged intermediates and (ii) reduction of reorganization energy. Although both of these effects cause a lowering of activation barriers, it is appropriate to distinguish them from each other because the ``mechanism'' by which the barriers are reduced are fundamentally different in the two cases. In the first case, the transition states flanking an intermediate are lowered as a consequence of the stabilization of that intermediate (Fig. 5B). This is basically the rationale behind so-called linear free energy relationships, such as the Hammond postulate, where one observes a systematic dependence of DeltaGon DeltaG^o for a series of similar reactions. The second effect, on the other hand, corresponds to a ``pure'' transition state stabilization that does not affect DeltaG^o for the given reaction step (Fig. 5C). The intervening barrier between two states is instead lowered as a consequence of a reduced energy gap between the potential surfaces of these states. In the terminology of electron transfer theory(30) , the reorganization energy () for a diabatic reaction is defined as the (absolute) free energy change that would be measured on the product potential surface if the geometry, or reaction coordinate, of the system is changed from the reactant surface minimum to the product minimum (Fig. 5). For adiabatic reactions, such as the proton transfers considered here, the reorganization energy can still be defined in terms of the diabatic surfaces, although the barrier height is not solely determined by the intersection of these but also by the magnitude of the off-diagonal matrix elements(31) . That is to say that the quantum mechanical coupling between the states causes a lowering of the activation barrier compared with the diabatic case where this coupling is zero.


Figure 5: Schematic representation of the two observed sources of catalysis. A, uncatalyzed reaction. B, catalysis by reaction energy reduction (e.g. stabilization of intermediates). C, catalysis by reorganization energy reduction. The two solid curves in each diagram are the diabatic reactant and product free energy surfaces, and the dashed curve represents the adiabatic ground state surface. The reorganization energy is denoted by .



The first of the above effects favors the transfer of negative charge from Glu to the substrate as well as His. The proximity of Lys (which also forms an ion pair with Glu) to the substrate together with hydrogen bonds from His and Asn are found to be important for the stability of the enediolate species, in agreement with earlier conjectures(6, 9) . Lodi and Knowles (32) have also pointed out the presence of the ``well aimed'' alpha-helix directed at His that presumably contributes to the stability of the imidazolate ion. Furthermore, the present simulations indicate that an active site water molecule, located approximately between Glu and His (denoted 626 in the PDB entry 7tim), plays an important role during catalysis (Fig. 6). This water has sufficient rotational and translational freedom so as to follow the ``negative charge'' because it is translocated within the active site, and it can provide stabilization by hydrogen bonding to Glu, the substrate, and His when these are negatively charged. In fact, the H-bonding network in the active site appears to be particularly well suited for adaptation to the movement of negative charge during catalysis.


Figure 6: Stereo snapshot of the TIM active site in the intermediate state ((3)) with His negatively charged. The water molecule corresponding to number 626 in the PDB entry 7tim can be seen between His and Glu. Other water molecules in the active site and dipolar NH groups of the protein solvating the phosphate anion are also shown.



Based on the above observation regarding the role of water one would expect that a mutation that displaces this water molecule should hamper catalysis. Judging from the structure mutations of Ser, Cys, or Ile might cause such an effect. Interestingly, we note that Blacklow and Knowles (33) have examined the S96P mutant and found that it considerably reduces the catalytic efficiency. This mutant would indeed appear to displace water, and our results suggest that this could be the reason for its reduced activity. Furthermore, we observe that the hydroxyl group of Thr of the second subunit can interact with the deprotonated histidine in the (3) intermediate state (Fig. 6). As in the crystal structure, Thr donates a hydrogen bond to the carboxylate group of Glu when His is neutral, but in the second reaction step the hydrogen bond can switch to the negative imidazolate ring of His, thereby providing additional stabilization of this group. This behavior could thus provide a new clue to the low enzymatic activity of monomeric TIM (34) , because Thr belongs to the second subunit. Borchert et al.(34) have shown by protein engineering of trypanosomal TIM that a shortening of the loop around Thr (involving deletion of this residue) yields a stable monomeric enzyme. This TIM variant still retains most of its substrate affinity (20-fold lower for LGAP than the wild type) while having a 1000-fold lower value of k(34) . Thr (from the second subunit) has also been implicated as an important residue in the work by Brown and Kollman(35) , for a somewhat different reason. These authors carried out a 5 ps MD simulation of monomeric TIM in vacuum and found considerable disruptions of the active site as a result.

Our simulations show that the mechanism proposed by Bash et al.(9) involving the imidazolate species is energetically and structurally very reasonable. The distinction between this mechanism and that proposed by Alagona et al.(10, 13) , which involves intramolecular proton transfer between the substrate oxygens, is, however, rather subtle. The difference only pertains to whether a proton is relayed via His or whether the -NH dipole on this residue merely provides a hydrogen bond during intramolecular transfer. The barrier obtained by Alagona et al.(13) in vacuum for this process is only about 14 kcal/mol, and more recent calculations including protein residues within 10 Å of the substrate and a realistic model of DHAP yielded a barrier of about 12 kcal/mol(10) . One can thus not exclude the possibility that the enzyme also could function via such an intramolecular process. A possible problem with this alternative, however, is that the His dipole would seem to interact less favorably with the intramolecular transition state than with the enolates. Comparing Fig. 11 (c and d) of (10) also suggests that this is the case, because the barrier increases by some 5 kcal/mol when His (plus all other residues within 10 Å of the substrate) are included as opposed to the calculation with only Asn, Lys, Ser, and Glu present(10) .

In the present work we find that the barriers for the intermolecular mechanism are low enough not to be rate determining. This is in contrast with (9) and (10) where they are too high to be compatible with experiment. It should also be noted here that these studies did not calculate room temperature free energies but rather minimum energy paths at zero temperature, which may involve more serious problems with local minima. Furthermore, in Ref.10 water molecules do not appear to have been included, which could lead to exaggerated electrostatic effects. In our simulations we also find that the pK difference between imidazole and the enediol is not altered very much by the enzymic environment. However, in both the work of Alagona et al.(10) and that of Bash et al.(9) , the proton transfer from imidazole is considerably ``uphill.'' This might seem a bit awkward in view of the experimental finding that the (first) pK of His is significantly lowered in enzyme(8) . Alagona et al.(10) further argue that the fact that the nearby lysine residue (Lys) normally has a lower pK than imidazole would disfavor His as a proton donor. This may, however, be of less relevance because (i) the pK of His is known to be lowered(8) , as mentioned above, and (ii) Lys participates in a salt bridge with Glu, which presumably raises its pKvalue. Regarding the issue of inter- versus intramolecular proton transfer pathways, it has also been shown that the mutant enzyme with His replaced by Gln, rather than to employ intramolecular proton transfer, uses Glu to effect all the proton tranfers(36) . However, we can of course not make any judgement of the energetics of the intramolecular process on the basis of the calculations reported here. On the other hand, our simulations suggest that the proton transfer barriers of the imidazole mechanism are at least low enough to not limit the rate of the enzyme.

The ability to reduce reorganization energies seems to be an important property of enzymes(15, 37) . For a given chemical reaction step, part of the energy barrier can be viewed as an intrinsic contribution that does not depend on DeltaG^o (the reaction free energy) but on the curvature of the free energy functions of the reactant and product states, which reflects the reorganization, e.g. of dipolar groups, involved in the reaction(30) . In the case of large energy gaps between the reactant and product surfaces at their respective minima, there will thus be a considerable reorganization of the system as the reaction proceeds, which gives rise to the intrinsic barrier. This is typically the case for charge separation and transfer reactions in aqueous solution, where water molecules must change their average polarization direction. In enzymes, polar groups are attached to a relatively rigid framework and are therefore less reorientable than in water. This fact, together with a clever design of the framework in which dipoles are pointing in a favorable direction, can lead to a considerable reduction of the reorganization energy compared with solution reactions. In the present calculations we observe this effect in all the steps of the catalyzed reaction. Fig. 7compares the diabatic free energy curves (i.e. the free energy curves of the pure VB states without considering their mixing) in water and in the enzyme for the states (2) and (3), corresponding to the second reaction step. These free energy curves are obtained by the same potential of mean force formulation as that used to calculate the ground state profile(31) . The reorganization energy upon going from (2) to (3) can be defined as = Delta - DeltaG^o, where Delta = (1) - (2) is the energy gap between the two curves at the minimum of (2). Vice versa, = Delta - DeltaG^o for the opposite reaction, Delta now being measured at the minimum of (3). In the case of the second reaction step, the absolute value of DeltaG^o is small (about 3 kcal/mol both in solution and enzyme), and the reduction of reorganization energy can therefore be seen directly from the fact that the energy gaps at the minima of the diabatic ``reactant'' and ``product'' free energy surfaces for the reaction step are clearly smaller in the enzyme than in solution. In this case, as can be seen from Fig. 4, which gives the actual adiabatic energies (i.e. the ground state free energy surface after mixing the VB states), the enzyme leaves DeltaG^o close to its uncatalyzed value and mainly affects the barrier height.


Figure 7: Diabatic free energy functions for the states (2) and (3) in water (upper curve) and in the enzyme (lower curve) as a function of the energy gap between the two states. The reorganization energies (see text for definition) upon transfer from the reactant to the product state and vice versa for this (the second) reaction step are indicated. The reduction of reorganization energy in the protein is reflected by the energy gaps being significantly smaller than in the uncatalyzed reaction.



Another convincing proof of the reorganization energy reduction effect can be obtained by artificially adjusting DeltaG^o (by varying the gas phase energy difference discussed above) for a given reaction step in water so that it coincides with the calculated value in the enzyme and then compare the resulting hypothetical barrier height with that in the enzyme. For the first step, one then obtains a barrier of 17.3 kcal/mol in water with the same reaction free energy as in the protein. The fact that the barrier obtained in this way is considerably higher than in the protein shows that the enzyme indeed does more than simply change relative pKs (the same conclusion is reached for the other steps as well). The stabilization of transition states can thus be achieved both by lowering the energies of intermediates and by reduction of the above-mentioned energy gaps or reorganization energies. One can, in principle, imagine several ways of achieving the latter, but at least when the reaction involves polar states it is clear that a protein matrix with properly oriented dipoles (or electric fields) can provide a microenvironment that responds with a smaller reorganization than water while still providing a substantial solvation of polar states. It is interesting to note here that the possibility of reorganization energy reduction as a source of catalysis in enzymes was conjectured already in 1980 by Albery(38) .

In view of the observations above it is interesting to try to examine Menger's recently proposed ``split-site model'' with which he has attempted to rationalize the catalytic effect of TIM without invoking TS stabilization(39) . Menger instead suggests that binding energy at a site (viz. the phosphate group) at some distance from the reaction centre is used to destabilize the reactant ground state by forcing the substrate into a conformation characterized by unfavorable ``compressive and desolvation forces'' that raise the reactant energy level and accelerate the reaction (see also the review by Page (40) and references therein). In our simulations we find no evidence of substrate compression or strain, the average intramolecular substrate energies being very similar in water and in the protein. As far as the desolvation idea is concerned, its implication here is a possible destabilization (upwards pK shift) of the negatively charged Glu upon substrate binding. There is no doubt that the pKdifference between DHAP and the catalytic base is reduced by the enzyme (cf. Fig. 4) and to some extent the question of where the absolute pKs end up is of less importance. We can, however, make the following observations. The negatively charged intermediates of the reaction are clearly stabilized by several hydrogen bonds, and their electrostatic interaction energy with the surroundings is considerably more favorable than in water. That the active site is biased toward negative charge is also evidenced by the observed (8) pK drop of His and by the strong binding of the phosphate group. A comparison of the evironment around Glu in the unliganded enzyme (5) and the substrate reactant state (R) shows an equal number of three hydrogen bonds to the carboxylate group (in the latter case from the thiol group of Cys, the DHAP hydroxyl group, and a water molecule). Thus, the most obvious explanation for the small reaction free energy differences of the first and last steps is that there is large stabilization of the negatively charged intermediates. This does not, of course, exclude the possibility of the Glu pK being shifted somewhat upwards but the main effect does not appear to be ``ground state destabilization.'' Instead, our results clearly implicate the two effects mentioned above as responsible for the catalytic power of TIM.

It should perhaps again be emphasized that the two general effects implicated here as the source of catalysis are not in any way at variance with the notion of transition state stabilization. Rather they constitute two different ways of achieving such a stabilization, and in that sense provides a clearer picture of what the enzyme does. It seems very likely to us that these two effects constitute general principles for enzymic catalysis of charge transfer and separation reactions. Certainly, more investigations are required to establish this point, but a considerable amount of data has accumulated during recent years that seem to support such a conjecture(14, 15, 17, 31, 37) .


FOOTNOTES

*
This work was supported by the Swedish Natural Science Research Council. The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore by hereby marked ``advertisement'' in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

§
To whom correspondence should be addressed. Tel.: 46-18174109; Fax: 46-18536971.

Recipient of a post-doctoral fellowship from the foundation Wenner-Grenska Samfundet.

(^1)
The abbreviations used are: TIM, triosephosphate isomerase; DHAP, dihydroxyacetone phosphate; GAP, glyceraldehyde 3-phosphate; MD, molecular dynamics; FEP, free energy perturbation; EVB, empirical valence bond; VB, valence bond.

(^2)
J. Åqvist, unpublished data.


ACKNOWLEDGEMENTS

We thank Prof. Donald Truhlar for kindly making the AMSOL program available to us and Prof. Greg Petsko for sending us the x-ray coordinates of the TIM-PGH complex. We also thank Tomas Hansson for useful discussions.


REFERENCES

  1. Rieder, S. V., and Rose, I. A. (1959) J. Biol. Chem. 234, 1007-1010 [Free Full Text]
  2. Knowles, J. R., and Albery, W. J. (1977) Acc. Chem. Res. 10, 105-111
  3. Nickbarg, E. B., and Knowles, J. R. (1988) Biochemistry 27, 5939-5947 [Medline] [Order article via Infotrieve]
  4. Knowles, J. R. (1991) Nature 350, 121-124 [CrossRef][Medline] [Order article via Infotrieve]
  5. Lolis, E., Alber, T., Davenport, R. C., Rose, D., Hartman, F. C., and Petsko, G. A. (1990) Biochemistry 29, 6609-6618 [Medline] [Order article via Infotrieve]
  6. Davenport, R. C., Bash, P. A., Seaton, B. A., Karplus, M., Petsko, G. A., and Ringe, D. (1991) Biochemistry 30, 5821-5826 [Medline] [Order article via Infotrieve]
  7. Richard, J. P. (1984) J. Am. Chem. Soc. 106, 4926-4936
  8. Lodi, P. J., and Knowles, J. R. (1991) Biochemistry 30, 6948-6956 [Medline] [Order article via Infotrieve]
  9. Bash, P. A., Field, M. J., Davenport, R. C., Petsko, G. A., Ringe, D., and Karplus, M. (1991) Biochemistry 30, 5826-5832 [Medline] [Order article via Infotrieve]
  10. Alagona, G., Ghio, C., and Kollman, P. A. (1995) J. Am. Chem. Soc. 117, 9855-9862
  11. Daggett, V., Brown, F., and Kollman, P. (1989) J. Am. Chem. Soc. 111, 8247-8256
  12. Joseph, D., Petsko, G. A., and Karplus, M. (1990) Science 249, 1425-1428 [Medline] [Order article via Infotrieve]
  13. Alagona, G., Desmeules, P., Ghio, C., and Kollman, P. A. (1984) J. Am. Chem. Soc. 106, 3623-3632
  14. Warshel, A. (1991) Computer Modeling of Chemical Reactions in Enzymes and in Solutions , John Wiley & Sons, Inc., New York
  15. Åqvist, J., and Warshel, A. (1993) Chem. Rev. 93, 2523-2544
  16. van Gunsteren, W. F., and Berendsen, H. J. C. (1987) Groningen Molecular Simulation (GROMOS) Library Manual , Biomos B.V., Nijenborgh, Groningen, The Netherlands
  17. Åqvist, J., Fothergill, M., and Warshel, A. (1993) J. Am. Chem. Soc. 115, 631-635
  18. Cramer, C. J., and Truhlar, D. G. (1992) Science 256, 213-217
  19. Åqvist, J. (1990) J. Phys. Chem. 94, 8021-8024
  20. Pauling, L. (1945) The Nature of the Chemical Bond , Cornell University Press, Ithaca, New York
  21. Chiang, Y., Kresge, A. J., Tang, Y. S., and Wirz, J. (1984) J. Am. Chem. Soc. 106, 460-462
  22. Guthrie, J. P. (1979) Can. J. Chem. 57, 1177-1185
  23. Eigen, M. (1964) Angew. Chem. Int. Ed. Engl. 3, 1-72
  24. Eigen, M., and Hammes, G. G. (1963) Adv. Enzymol. 25, 1-38 [Medline] [Order article via Infotrieve]
  25. Jorgensen, W. L. (1989) Acc. Chem. Res. 22, 184-189
  26. Beveridge, D. L., and DiCapua, F. M. (1989) Annu. Rev. Biophys. Biophys. Chem. 18, 431-492 [CrossRef][Medline] [Order article via Infotrieve]
  27. Straatsma, T. P., and McCammon, J. A. (1992) Annu. Rev. Phys. Chem. 43, 407-435 [CrossRef]
  28. Åqvist, J., and Warshel, A. (1992) J. Mol. Biol. 224, 7-14 [Medline] [Order article via Infotrieve]
  29. King, G., and Warshel, A. (1989) J. Chem. Phys. 91, 3647-3661 [CrossRef]
  30. Cohen, A. O., and Marcus, R. A. (1968) J. Phys. Chem. 72, 4249-4256
  31. Warshel, A., Hwang, J.-K., and Åqvist, J. (1992) Faraday Discuss. Chem. Soc. 93, 225-238
  32. Lodi, P. J., and Knowles, J. R. (1993) Biochemistry 32, 4338-4343 [Medline] [Order article via Infotrieve]
  33. Blacklow, S. C., and Knowles, J. R. (1990) Biochemistry 29, 4099-4108 [Medline] [Order article via Infotrieve]
  34. Borchert, T. V., Abagyan, R., Jaenicke, R., and Wierenga, R. K. (1994) Proc. Natl. Acad. Sci. U. S. A. 91, 1515-1518 [Abstract]
  35. Brown, F. K., and Kollman, P. A. (1987) J. Mol. Biol. 198, 533-546 [Medline] [Order article via Infotrieve]
  36. Nickbarg, E. B., Davenport, R. C., Petsko, G. A., and Knowles, J. R. (1988) Biochemistry 27, 5948-5960 [Medline] [Order article via Infotrieve]
  37. Yadav, A., Jackson, R. M., Holbrook, J. J., and Warshel, A. (1991) J. Am. Chem. Soc. 113, 4800-4805
  38. Albery, J. W. (1980) Annu. Rev. Phys. Chem. 31, 227-263
  39. Menger, F. M. (1992) Biochemistry 31, 5368-5373 [Medline] [Order article via Infotrieve]
  40. Page, M. I. (1987) in Enzyme Mechanism (Page, M. I., and Williams, A., eds), Chapter 1, Royal Society of Chemistry, London

©1996 by The American Society for Biochemistry and Molecular Biology, Inc.