Computer Simulation of the Triosephosphate Isomerase Catalyzed Reaction*

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 triosephos-phate 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 re-organization energy, which are likely to be general principles of enzyme catalyzed charge transfer and separation reactions.

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.
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 165 , yielding a corresponding enediolate species. (ii) This enediolate is then protonated at O-2 by the neutral imidazole ring of His 95 , giving a doubly protonated enediol and an imidazolate anion. (iii) His 95 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 165 , 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 Ϫ1 (3,7). As a comparison the rate-limiting step for nonenzymatic conversion of DHAP is about 10 5 s Ϫ1 , in which case the mechanism involves intramolecular proton abstraction by the phosphate group (7).
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 2Ϫ moiety was, however, represented with the standard harmonic bond functions (16) because these bonds are not broken during the reaction.
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, ␣, 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 me-chanical 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 I), 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 (␣ values) between the VB states and the off-diagonal matrix elements H ij (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 a values of hydrogens on C-1 and C-2 as well as the associated activation barriers. Starting from the pK a of acetone in water of 19.2 (21), these pK a s can, however, be estimated by the substituent effects of a hydroxyl group and the phos-2 J. Åqvist, unpublished data.  (16) were used except for the oxygens bearing a large negative charge (͉q͉ Ͼ 0.5) for which the parameters derived in Ref. 17 were used (these have been calibrated to reproduce hydration free energies for carboxylate ions). phate 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 a 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 a 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 a 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 a at 17.2 is also rather close to the value inferred by Richard (7)). The corresponding pK a value of 15.9 for LGAP is related by the observed equilibrium constant of 1 ⁄22 to that of DHAP (7). By interpolation in the diagram, the free energy relationship of Fig. 3 can also be used to deduce the activation barriers for proton exchanges with the carboxylate moiety of glutamic acid in aqueous solution (pK a ϭ 4.1), which are used as reference reactions for the first and fourth step of the enzymic conversion (cf. Table I). 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 95 . The pK a difference between donor and acceptor in solution is about 2 pK units, with the neutral imidazole ionizing at the higher pK a 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 Ref. 23 and Table V of Ref. 24), and the corresponding energetics are listed in Table I.
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, ␣, between the VB states and of the off-diagonal hamiltonian matrix elements (the latter are represented by functions of the form ij ϭ A ij e Ϫ ij rxy , where r XY denotes the distance between the two heavy atoms that exchange the proton (14, 15)). Table I lists the EVB parameters obtained from this calibration procedure together with the corresponding reaction free energies and barriers that were estimated for the solution reactions. After performing the initial calibration simulations in water the obtained parameters are kept fixed and the corresponding reaction steps are simulated inside the enzyme. Here, it should be pointed out that we neglect possible differences in barrier heights and reaction free energies between the solution and enzyme reactions that are due to quantum mechanical effects such as proton tunneling and zero-point energies.
That is, the calibration in water yields an effective potential surface that includes such effects, but to the extent that quantum corrections to the classical free energies differ between the enzyme and solution reactions they are neglected. Furthermore, we only examine the nonconcerted mechanism here and assume that off-diagonal matrix elements H ij for which ͉i Ϫ j͉ Ն 2 (indicative of simultaneous proton tranfers) are zero. The simulation characteristics were the same as in Refs. 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 Ref. 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. 3 and 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  Ref. 22, squares) and on base catalyzed phosphate elimination reactions (7) of DHAP (circles) and LGAP (triangles) using the estimated pK a values of DHAP and LGAP. The latter pK a was related to the former by the observed equilibrium constant of 1 ⁄22 between LGAP and DHAP in solution (7). All values have been corrected for the number of equivalent protons.

TABLE I EVB parameters and estimated energetics for the uncatalyzed reference reactions in water used for calibration
All quantities are in kcal/mol except ij , which is given in Å Ϫ1 . ⌬G aq o and ⌬G aq ‡ measure the reaction and activation energies, respectively, for each reaction step in the forward direction (cf. Fig. 1). Note that the gas phase energy differences (see text), ␣, cannot directly be equated to the true free energy difference between the reactant and product fragments of each step (as in Ref. 17), because they in this case implicitly include nonbonded interactions within each fragment that may be non-zero in the equilibrium conformation.
Step  (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 ⌬G ‡ on ⌬G 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 ⌬G 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.
The first of the above effects favors the transfer of negative charge from Glu 165 to the substrate as well as His 95 . The proximity of Lys 12 (which also forms an ion pair with Glu 97 ) to the substrate together with hydrogen bonds from His 95 and Asn 10 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" ␣-helix directed at His 95 that presumably contributes to the stability of the imidazolate ion. Furthermore, the present simulations indicate that an active site water molecule,  Fig. 4

assuming that only the lowest energy intermediate is kinetically observable
Step  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. located approximately between Glu 165 and His 95 (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 165 , the substrate, and His 95 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.
Based on the above observation regarding the role of water 626 one would expect that a mutation that displaces this water molecule should hamper catalysis. Judging from the structure mutations of Ser 96 , Cys 126 , or Ile 170 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 626 , and our results suggest that this could be the reason for its reduced activity. Furthermore, we observe that the hydroxyl group of Thr 75 of the second subunit can interact with the deprotonated histidine in the ⌿ 3 intermediate state (Fig. 6). As in the crystal structure, Thr 75 donates a hydrogen bond to the carboxylate group of Glu 97 when His 95 is neutral, but in the second reaction step the hydrogen bond can switch to the negative imidazolate ring of His 95 , 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 75 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 75 (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 cat (34). Thr 75 (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 95 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 95 dipole would seem to interact less favorably with the intramolecular transition state than with the enolates. Comparing Fig. 11 (c and d) of Ref. 10 also suggests that this is the case, because the barrier increases by some 5 kcal/mol when His 95 (plus all other residues within 10 Å of the substrate) are included as opposed to the calculation with only Asn 10 , Lys 12 , Ser 96 , and Glu 97 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 Refs. 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 a 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 a of His 95 is significantly lowered in enzyme (8). Alagona et al. (10) further argue that the fact that the nearby lysine residue (Lys 12 ) normally has a lower pK a than imidazole would disfavor His 95 as a proton donor. This may, however, be of less relevance because (i) the pK a of His 95 is known to be lowered (8), as mentioned above, and (ii) Lys 12 participates in a salt bridge with Glu 97 , which presumably raises its pK a value. Regarding the issue of inter-versus intramolecular proton transfer pathways, it has also been shown that the mutant enzyme with His 95 replaced by Gln, rather than to employ intramolecular proton transfer, uses Glu 165 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 ⌬G 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. 7 compares 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 233 ϭ ͉⌬⑀͉ Ϫ ⌬G 233 o , where ⌬⑀ ϭ ⑀ 1 Ϫ ⑀ 2 is the energy gap between the two curves at the minimum of ⌿ 2 . Vice versa, 332 ϭ ͉⌬⑀͉ Ϫ ⌬G 332 o for the opposite reaction, ⌬⑀ now being measured at the minimum of ⌿ 3 . In the case of the second reaction step, the absolute value of ⌬G 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 ⌬G o close to its uncatalyzed value and mainly affects the barrier height.
Another convincing proof of the reorganization energy reduction effect can be obtained by artificially adjusting ⌬G 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 FIG. 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. than in the protein shows that the enzyme indeed does more than simply change relative pK a s (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 a shift) of the negatively charged Glu 165 upon substrate binding. There is no doubt that the pK a difference between DHAP and the catalytic base is reduced by the enzyme (cf. Fig. 4) and to some extent the question of where the absolute pK a s 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 a drop of His 95 and by the strong binding of the phosphate group. A comparison of the evironment around Glu 165 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 126 , 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 165 pK a 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 stabiliza-tion. 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).