HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Full Text (PDF) Alert me when this article is cited Alert me if a correction is posted Citation Map Services Email this article to a friend Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Garcia-Viloca, M. Articles by Gao, J. Search for Related Content PubMed PubMed Citation Articles by Garcia-Viloca, M. Articles by Gao, J. Social Bookmarking                   What's this?

Sensitivity of molecular dynamics simulations to the choice of the X-ray structure used to model an enzymatic reaction

Department of Chemistry and Supercomputer Institute, University of Minnesota, Minneapolis, Minnesota 55455, USA

Reprint requests to: Donald G. Truhlar, Department of Chemistry and Supercomputer Institute, University of Minnesota, Minneapolis, MN 55455, USA; e-mail: truhlar{at}umn.edu; fax: (612) 626-9390.

(RECEIVED December 13, 2003; FINAL REVISION January 19, 2004; ACCEPTED January 20, 2004)

	Abstract

TOP Abstract Introduction Materials and methods Results and Discussion References

 
A subject of great practical importance that has not received much attention is the question of the sensitivity of molecular dynamics simulations to the initial X-ray structure used to set up the calculation. We have found two cases in which seemingly similar structures lead to quite different results, and in this article we present a detailed analysis of these cases. The first case is acyl-CoA dehydrogenase, and the chief difference of the two structures is attributed to a slight shift in a backbone carbonyl that causes a key residue (the proton-abstracting base) to be in a bad conformation for reaction. The second case is xylose isomerase, and the chief difference of the two structures appears to be the ligand sphere of a Mg²⁺ metal cofactor that plays an active role in catalysis.

Keywords: combined quantum mechanics/molecular mechanics; molecular dynamics; potential of mean force; structure-based enzyme modeling

Article and publication date are at http://www.proteinscience.org/cgi/doi/10.1110/ps.03504104.

	Introduction

TOP Abstract Introduction Materials and methods Results and Discussion References

 
Molecular dynamics (MD) simulations of enzyme reactions begin with the construction of a model for the Michaelis complex, typically making use of X-ray crystal structures that are complexed with a substrate analog. Except in a few rare cases, the X-ray structures are not the Michaelis complex with the native substrate bound, which would be the most relevant structure for the simulation. Sometimes the available structures are obtained from a different species (e.g., bacterial rather than human), or the enzyme has no ligand or is bound to an inhibitor or product, or the enzyme has been mutated or modified to facilitate the crystallization process. In other cases the crystal may have different quaternary structure (e.g., monomer, dimer, tetramer) than the active form of the enzyme. Any of these variations may significantly or subtly alter the structure and function of the native biological system. Furthermore, in most cases the experimental conditions (such as pH, temperature, and co-solute) used to crystallize proteins are very different from the biological conditions under which a given enzyme functions in vivo. Little is known of the effect of these experimental conditions on the enzyme structure, and it is possible that attendant changes in the enzyme could mask critical interactions necessary for its in vivo efficiency. Often there are two or more structures of a given enzyme that differ in binding partner or in crystallization conditions. In this article we consider a question that is relevant to all structure-based computer studies of enzymatic reactions, in particular, the importance of the choice of X-ray structure that is used to obtain initial coordinates for the simulation.

In the present article we present two examples, acyl-CoA dehydrogenase (ACAD)¹ and xylose isomerase (XyI),² in which the initial X-ray structure chosen to model an enzymatic reaction was found to be in a significantly less reactive conformation than a second choice made later. In particular, only the use of the second structure allowed us to build a model that yielded reasonable results for the enzymatic reaction studied. We hope that these examples illustrate a crucial element in the study of enzymatic reactions, even before the simulation has started, namely, the criticality of the choice of X-ray structure that is used to build a model for the Michaelis complex. For each enzyme we will compare potentials of mean force (PMFs) computed using two different X-ray structures as starting point.

	Materials and methods

TOP Abstract Introduction Materials and methods Results and Discussion References

 
We calculate PMFs (also called free energy profiles) by using umbrella sampling techniques in MD simulations based on combined quantum mechanical–molecular mechanical (QM/MM; Warshel and Levitt 1976; Singh and Kollman 1986; Bash et al. 1987; Gao 1995; Gao et al. 1998; Gao and Truhlar 2002) potential energy functions. The computational details have been given in previous articles on these two enzymatic systems (Garcia-Viloca et al. 2002, 2003; Poulsen et al. 2003). In both studies, the active site is represented quantum mechanically with a semiempirical method (Dewar et al. 1985; Stewart 1989). For the electrostatics, we include all charge–charge interactions within a cutoff.

The original aim of these studies was to elucidate the catalytic mechanism of the ACAD and XyI enzymes and to calculate the primary kinetic isotope effect on the rate constants. For the latter objective, we needed to improve the accuracy of the QM/MM potential energy surfaces, and we used (Garcia-Viloca et al. 2003; Poulsen et al. 2003) specific reaction parameters (Gonzalez-Lafont et al. 1991) for the semi-empirical method or we added a valence bond correction term (Alhambra et al. 1999; Devi-Kesavan et al. 2003) to the total potential energy. In the present article, however, for the purpose of comparison, all calculations were carried out without specific reaction parameters or valence bond terms. Therefore, the free energy profiles are not quantitatively identical to those of Garcia-Viloca et al. (2003) and Poulsen et al. (2003).

In the simulations of ACAD, 55 atoms were treated quantum mechanically by the semiempirical Austin model 1 (AM1; Dewar et al. 1985). These atoms include (Poulsen et al. 2003) part of the substrate, the flavin adenine dinucleotide (FAD) cofactor, and the side-chain of Glu376 for medium-chain ACAD or the side-chain of Glu367 for butyryl-CoA dehydrogenase (BCAD). The rest of the protein-solvent system was represented classically by the CHARMM (MacKerell et al 1998) force field and the TIP3P (Jorgensen et al. 1983) water model. For the proton and the hydride transfer steps, the reaction coordinate was defined as the difference between the distance of the proton or hydride atom from the donor atom and its distance from the acceptor atom (Poulsen et al. 2003).

In the simulations of the XyI reaction, all 19 atoms of the substrate (D-xylose) were represented quantum mechanically by the semiempirical PM3 method (Stewart 1989). The enzymatic residues and the solvent were represented by the CHARMM force field and the TIP3P water model, respectively. The reaction coordinate was defined as the difference between the breaking bond distance and the forming bond distance of the hydride transfer step (Garcia-Viloca et al. 2003).

In all simulations in this article, we use stochastic boundary conditions (Brooks et al. 1985) in which the atoms in a spherical shell centered on the active site are treated by restrained Langevin dynamics. The shell extends from 25 Å to 30 Å for ACAD and from 20 Å to 24 Å for XyI; the region inside the inner radius is called the reaction zone and is treated by unrestrained Newtonian dynamics. All atoms in the stochastic buffer shell are restrained by a harmonic restoring potential to remain close to their original positions, and in addition, they are subject to Langevin forces, which consist of a frictional force and a thermal random force. A consequence of the restraint is that certain large-amplitude protein motions are excluded. For example, if there is a hinge or domain motion in the tertiary or quaternary structure, the hinge might lie outside of the reaction zone. Thus, the calculations involve the implicit assumption that large-amplitude motions of this extent are not important for the chemical step (e.g., we assume that they do not couple to the reaction coordinate), although, for example, they might be important for substrate binding or product release. A more subtle possibility is that such large-amplitude motion might be coupled to conformational transitions near the active site, and excluding them could make it more difficult for the simulation to evolve into the most productive conformation. Local conformational changes of a few tenths to 2 Å that can be accomplished with small barriers are taken into account, and this may be very important.

	Results and Discussion

TOP Abstract Introduction Materials and methods Results and Discussion References

 
Acyl-CoA dehydrogenase
The family of ACADs is very large, with nine members (Ghisla and Thorpe 2004). We are concerned here with the subfamily that catalyzes the first step of the -oxidation cycle of straight-chain fatty acids; this step involves the oxidation of fatty acyl-CoA substrates to trans-enoyl-CoA products.

Various members of subclass of straight-chain dehydrogenating enzymes may be considered, and they have different substrate specificities. We consider short (SCAD) and medium (MCAD) chain ACADs. Although the chain-length specificities depend on pH and other conditions, SCADs are typically optimal for fatty acids with chain lengths of about four carbon atoms, and MCADs are typically optimal for chain lengths of about eight (Ikeda et al. 1985; Finocchario et al. 1987; Trivel et al. 1995; Nandy et al. 1996; Ghisla et al. 2003), but the specificity ranges tend to be fairly broad and overlapping (e.g., both SCADs and MCADs are active for chain lengths in the range of four to six). In the present study, we compare two MD simulations, one based on the X-ray structure of mammalian MCAD (Lee et al. 1996) and the other based on BCAD (Djordjevic et al. 1995), which is a bacterial analog of the SCAD from mammalian mitochondria. BCAD can be purified in high yield from bacterial cells, and therefore, a number of experimental studies have been performed on BCAD. In Figure 1, we show the FAD cofactor, the substrate, and the participating enzymatic residue along with a proposed mechanism. A body of evidence summarized elsewhere (Ghisla et al. 1984; Powell and Thorpe 1988; Bross et al. 1990; Vock et al. 1998; Ghisla and Thorpe 2004) indicates that, in ACADs, a glutamate (Glu367 in BCAD and Glu376 in MCAD) initiates the oxidation process by abstracting the substrate -proton. The overall reaction involves the concomitant or subsequent transfer of a hydride ion from the -carbon to the N5 position of flavin, resulting in the reduction of the FAD. These steps yield the two-electron–reduced flavin (FADH₂ or FADH^–, which is the deprotonated form of FADH₂) and oxidized enoyl-CoA product (Palfrey and Massey 1998; Ghisla and Thorpe 2004). Conceivably, the reduction of FAD can proceed via two extreme mechanisms, by a concerted proton–hydride transfer or a stepwise process (which has two possibilities: proton, then hydride; or proton, then electron, then hydrogen atom). In the present article, we only consider the proton followed by hydride stepwise mechanism, and we focus on the hydride step. We refer to recent articles (Rudik et al. 1998; Vock et al. 1998; Engst et al. 1999; Tamaoki et al. 1999; Pellett et al. 2000; Peterson et al. 2001; Bach et al. 2002; Gopalan and Srivastava 2002; Lamm et al. 2002; Ghisla and Thorpe 2004; Poulsen et al. 2003) for additional discussion of the reaction mechanism.

View larger version (14K): [in this window] [in a new window]   Figure 1. Schematic representation of the oxidation of butyryl CoA. The figure shows two possible resonance structures (mesormeric forms) of the reduced flavin.

In the simulations of BCAD, we used the crystal structure with the PDB code 1BUC [PDB] (Djordjevic et al. 1995). This structure was given as a dimer, from which the tetrameric structure is obtained by symmetry transformation based on the space group of the crystal. The structure contains an inhibitor, acetoacetyl-CoA. Compared with the normal substrate, the inhibitor is blocked in the -position with a C = O group instead of the corresponding CH2 group in the substrate. We replaced the C = O group with a CH2 group to obtain a model for the Michaelis complex structure. To accommodate the change in the geometry of the substrate after the replacement, the geometry of the part of the system included in a sphere of 20 Å centered at the active site was optimized for 40 steps of energy minimization with the rest of the coordinates frozen at the X-ray geometry. Then stochastic boundary MD simulations were carried out (see Materials and Methods) in order to allow adjustment of the environment to the ligand change.

Experimental data for kinetic isotope effects are available for MCAD (Pohl et al. 1986; Schopfer et al. 1988). It is not clear a priori whether these can be or should be reproducible by simulations that do not take account of small differences in protein structures, and one goal of the present article is to present some of our findings relevant to this question.

Simulations using structures derived from MCAD and BCAD with n-octanoyl and butyryl substrates, respectively, will be called simulation 1 and simulation 2, respectively. In all simulations, we included the dynamics of all subunits of the tetramer that are within 30 Å of one of the active sites. This active site was used to model the chemical reaction, and this active site was solvated by a 30 Å sphere of previously equilibrated water molecules around the active center.

Structural differences in ACADs
Figure 2 shows a superimposition of the wild-type apoenzyme and the double-mutant protein, which contains the substrate n-octanoyl-CoA. The figure shows that the conformation of Glu99 is different between the wild type and the mutant. In the later structure, Glu255 and Glu99 are in close contact and appear to form a hydrogen bond. If this is the case, Glu255 and Glu99 must share a proton to maintain the short distance observed in the X-ray structure. In contrast, the wild-type enzyme does not contain a hydrogen bond between Thr255 and Glu99 in the apoenzyme.

View larger version (33K): [in this window] [in a new window]   Figure 2. Superimposition of the X-ray structures of the human wild-type protein (PDB: 1EGE [PDB] ; Lee et al. 1996) and the Glu376Gly/Thr255Glu double mutant (PDB: 1EGC [PDB] ; Lee et al. 1996). The former structure is shown in violet, whereas the color code for the mutant structure is as follows: green for carbon, red for oxygen, yellow for sulfur, and blue for nitrogen. The red labels indicate the mutated residues. The interaction distances of important hydrogen bonds are shown in angstroms. The arrows point out the carbon of enzymatic residues.

Simulations
The possible role of an equilibrium of conformational states of the activating glutamate has been discussed in previous work (Tamaoki 1999). The active site of the initial structure of MCAD from our first simulation is shown in Figure 3. The relative orientation of the basic residue Glu376 in the active site indicates that it is in a good position for reaction. The positively charged arginine (residue 256) is stabilized by hydrogen bonds to the backbone C = O group of Glu376 and to the side chain of Asp253. The system was heated and equilibrated for a total of 80 psec by running stochastic boundary MD at the temperature (277 K) used experimentally (Schopfer et al. 1988) to determine the rate constant and kinetic isotope effects. This provided the starting structure in simulation 1 for calculating the PMF by the umbrella sampling technique. During the equilibration stage, the carboxylate group of Glu376 was found to turn toward the charged Arg256 residue. This results in a strong interaction between Arg256 and the carboxylate group of Glu376, making Glu376 a much weaker base for proton abstraction. Concomitantly, both Asp253 and the backbone C = O of Glu376 moved away from Arg256. The resulting structure is depicted (as an instantaneous snapshot after equilibration, not a minimized structure) in Figure 4, which illustrates a strong electrostatic interaction between the ion pair of Glu376 and Arg256. This interaction hinders proton abstraction from the substrate and results in a very high barrier for the proton transfer. The structural change also affects the hydride transfer. Initially, the aliphatic chain of the substrate is sandwiched between the flavin ring and Glu376 with the -hydrogen oriented toward the carboxylate group of Glu376 and the -hydrogen in close proximity to the N5 atom of flavin. The formation of the ion pair between Glu376 and Arg256 in the distorted enzyme configuration pulls Glu376 away from this favorable conformation for the chemical step, and the subsequent proton transfer reaction necessarily forces the substrate closer to the new Glu376 position. Consequently, we found that the entire octanoyl-CoA migrates ~1 Å, placing the -hydrogen in a poor orientation for the hydride ion transfer.

View larger version (23K): [in this window] [in a new window]   Figure 3. The initial structure and relative orientation of the active center of simulation 1. The substrate is n-octanoyl-CoA, and key distances are given in angstroms.

View larger version (19K): [in this window] [in a new window]   Figure 4. The structure and relative orientation of the active center of simulation 1 after equilibration. The substrate is n-octanoyl-CoA, and key distances are given in angstroms.

View larger version (22K): [in this window] [in a new window]   Figure 5. The structure and relative orientation of the active center of simulation 2 after equilibration. The substrate is butyryl-CoA, and the unit of the shown distances is angstroms.

View larger version (12K): [in this window] [in a new window]   Figure 6. Classical potentials of mean force (PMFs) for the proton transfer obtained for simulation 1 (solid line) and simulation 2 (dotted line). By definition, each PMF is zero at the reactant state.

The strong interaction between Glu376 and Arg256 in simulation 1 also affects the PMF for the hydride transfer because the substrate moves away from the FAD ring, resulting in a poor conformation for hydride transfer. For this reaction step, the difference in the two reaction profiles (Fig. 7) is even greater than the proton transfer step. The barrier height for the hydride transfer in simulation 1 is 17 kcal/mole larger than the one in simulation 2.

View larger version (13K): [in this window] [in a new window]   Figure 7. Classical potentials of mean force for the hydride transfer obtained for simulation 1 (solid line) and simulation 2 (dotted line). Note that Figures 6 and 7 pertain to the same process.

The X-ray structures of MCAD show only small changes upon substrate binding (Kim et al. 1993). Four water molecules are displaced from the active site with the formation of two hydrogen bonds by the substrate thioester carbonyl with 2'-OH of the ribityl chain of the flavin, and with the backbone amide N-H of Glu376. In addition, the residues Glu99, Glu376, and Tyr375 move slightly. Glu99 turns by ~90°, making the substrate binding pocket deeper, and the OE1 and OE2 atoms of the Glu376 side chain move 1.1 Å and 3.78 Å, respectively (Kim et al. 1993).

The carboxyl group of Glu376 is intimately involved in modulating the microscopic environment of the active site of the enzyme during the course of ligand binding and catalysis (Bross et al. 1990; Nandy et al. 1996; Rudik et al. 1998; Peterson et al. 2001; Gopalan and Srivastava 2002).

Nandy et al. (1996) studied the kinetics of the double mutants Glu376Thr/Thr255Glu and Glu376Gly/Thr255Glu of MCAD, in which the abstracting base is moved to the position that it occupies in long-chain ACAD, as well as the Thr255Glu mutant, which differs from the wild type in charge and steric effects. They concluded that the proton abstraction has peculiar geometric requirements. They discussed the role of Arg256 but concluded that it is too far away to hydrogen bond to Glu376, which was confirmed by X-ray structures (Lee et al. 1996). Gopalan and Srivastava (2002) carried out experimental and molecular modeling studies on the Glu376Gln mutation that are relevant to the functions of the active site residues. They hypothesized that Glu376 is involved in structuring the cavity of the active site and the solvent and in modulating the dynamics of the protein conformational changes due to electrostatic interactions among Glu376, FAD, and the substrate. They surmised that electrostatic interactions between the negative carboxyl group of Glu376 and the partial positive charged isoalloxazine ring of the FAD maintain the active site of the enzyme and that the Glu376Gln substitution diminishes these interactions so that the active site cavity becomes more unstructured, resulting in a weaker protein–ligand interaction. The same argument applies to a Glu376Gly substitution.

The question of the of the other possible roles (e.g., electrostatic and steric) of Glu376, in addition to its function as the abstracting base, remains a burning question, as does the possible role of Arg256. Our results presented here throw fuel, not extinguisher, on this fire, and suggest that further studies of the roles of these residues would be valuable. They also provide a caution to modelers that subtle differences in starting structures of enzyme simulations might have large effects on the results. Clearly, we must be especially cautious in interpreting simulations that do not provide as much catalytic enhancement as experiment; this can apparently result from small defects in the structures used for the simulation and need not indicate that the postulated mechanism or the electronic structure description of the bond rearrangement itself is qualitatively incorrect.

One possible criterion for whether the enzyme structure should be considered reasonable is that, even though one may not perform full simulations on every step of the catalytic process, the enzyme structure has an arrangement of amino acids, ligands, and possible protonation states and conformational changes that allows one to propose a complete sequence of steps for the reaction mechanism. For example, we next turn to XyI, in which, although we modeled only one step in full, the structure allows a reasonable description of all the steps except the pyranose ring opening and closing steps, whereas another structure, with a different number of water molecules in the active site and different metal binding patterns, does not.

XyI: The use of an inhibitor-enzyme complex proposed to be a transition state analog
XyI catalyzes the interconversion of D-xylose and D-xylulose. In addition, the enzyme has the ability to isomerize D-glucose to D-fructose and for that reason has been extensively used in industry (Bhosale et al. 1996). The overall reaction consists of at least three chemical steps: (1) ring opening, (2) hydrogen transfer between C2 and C1, and (3) ring closure. The rate limiting step of the overall reaction is the intramolecular hydride shift from C2 to C1 (Collyer et al. 1990; Whitlow et al. 1991; Rangarajan and Hartley 1992; Zheng et al. 1993), which is accompanied by prior and post proton transfer from the O2 atom and to the O1 atom, respectively (Fig. 8). There are two magnesium ions in the active site of XyI that are required to achieve the enzyme function (Bhosale et al. 1996). The intramolecular hydride shift is relatively slow (kcat 0.6 to 18 sec^–1; Jenkins et al. 1992; Lambier et al. 1992). Several crystal structures of the substrate–enzyme complex, which have been determined in different laboratories (Farber et al. 1989; Collyer et al. 1990; Whitlow et al. 1991; Jenkins et al. 1992; Lavie et al. 1994), contain a mixture of the reactant and the product substrates. In the presence of the substrate (D-xylose or D-glucose), the electron density of one of the metals does not have a spherical shape, which does not allow the resolution of a unique position for the metal. As a result, in some X-ray structures of substrate–enzyme complexes, there are two coordinates for one of the magnesium ions, and it is generally accepted that this metal moves during the hydride-shift (Collyer et al. 1990; Whitlow et al. 1991; Jenkins et al. 1992; Lavie et al. 1994). On the other hand, a unique position of the metal has been determined in structures in which the enzyme was co-crystallized with a nonsubstrate molecule, in particular two structures resolved by Collyer et al. (1990) that were proposed to be analogs of the Michaelis complexes, and the structure determined by Allen et al. (1995) with an inhibitor bounded in the active site of XyI from Streptomyces olivochromogenes (PDB: 2GYI [PDB] ). The inhibitor, D-threonohydroxamic acid (THA), was designed to mimic the transition state of the isomerization step (Fig. 8) catalyzed by XyI. Lavie et al. (1994) have resolved the structure of two substrates, D-glucose and 3-O-methylglucose, ligated to the same XyI (PDB: 1XYB [PDB] and PDB: 1XYC [PDB] , respectively), and the structure of the enzyme without substrate (PDB: 1XYA [PDB] ). The tautomer of THA bound in the active center of XyI was deduced by Allen et al. (1995) from their structural data, and it is represented in Figure 9. These researchers showed that the strong binding of the inhibitor did not induce any gross conformational change, although the reported C RMS deviation of 0.27 Å for the enzyme main chain compared with the apoenzyme may be deceptively small (Lavie et al. 1994). On the basis of its high affinity and the structural similarities with the glucose complexed structure, Allen et al. (1995) postulated that THA resembles the proposed transition state for the enzyme-catalyzed hydride transfer reaction.

View larger version (8K): [in this window] [in a new window]   Figure 8. Schematic representation of the interconversion between D-xylose and D-xylulose.

View larger version (8K): [in this window] [in a new window]   Figure 9. Schematic representation of the transition state (TS) structure for the intramolecular hydride transfer in D-xylose molecule, and schematic representation of the inhibitor D-threonohydroxamic acid (THA).

View larger version (25K): [in this window] [in a new window]   Figure 10. Partial view of the active site of the X-ray structure with D-glucose bound to xylose isomerase (PDB: 1XYB [PDB] ; Lavie et al. 1994).

View larger version (26K): [in this window] [in a new window]   Figure 11. Partial view of the active center of the X-ray structure with the inhibitor THA bound to xylose isomerase (PDB: 2GYI [PDB] ; Allen et al. 1995).

On the basis of all the structural data provided by the four X-ray structures, Petsko and coworkers (Lavie et al. 1994; Allen et al. 1995) proposed the mechanism represented in Figure 12 for the proton and hydride transfer steps that follow after the ring opening of the substrate. In this mechanism, the Mg2 ligand OH1700 is the base that abstracts the C2 hydroxyl proton. Concomitant with this proton transfer or after it, Mg2 and some of its ligands move toward Mg1 due to increased attraction from the O2-alkoxide, arriving at the conformation represented by structure II in Figure 12. The migration of the Mg2 accompanying the proton abstraction is also associated with major changes in the ligand–metal interactions. In particular, Asp254 is substituted by a water molecule and the O2-alkoxide ion of xylose, and Asp254 is replaced byAsp256. The conformation represented by II is the proposed reactive conformation for the hydride transfer step (Allen et al. 1995). In contrast, Collyer et al. (1990) and Whitlow et al. (1991) had proposed a somewhat different mechanism. In their proposal, a water molecule ligated to the mobile metal becomes a hydroxide ion after transferring a proton to Asp256, which then becomes the base responsible for the proton transfer step. In agreement with Petsko’s mechanism (Lavie et al. 1994; Allen et al. 1995), only after the proton transfer step does the distance between the metals decrease, and the substrate forms two ligands to Mg2 (or Mn2 in the proposal of Whitlow et al. 1991), giving rise to the reactive conformation for the hydride transfer step (structure III in Fig. 12). However, in the Whitlow mechanism (Collyer et al. 1990; Whitlow et al. 1991) only the neutral Asp256 loses ligation to Mg2, which is replaced by the carbonyl group of the substrate, whereas Asp254 maintains one of its ligands to Mg2.

View larger version (27K): [in this window] [in a new window]   Figure 12. Schematic representation of the mechanisms proposed by Petsko and coworkers (IIIIV; Lavie et al. 1994; Allen et al. 1995) and Whitlow et al. (1991; IIIIV) for the proton and hydride transfer steps.

Having examined the X-ray structural information explained above, we chose the most recent X-ray structure that contains a transition state analog (XyI–THA complex, PDB: 2GYI [PDB] ) in the active site to build the Michaelis model. In fact, this is the only crystal structure in which the residues in the active site form coordination spheres to the Mg2 ion in the proposed reactive conformation for the hydride transfer reaction.

The XyI–THA structure contains the Cartesian coordinates for the enzyme atoms of the dimer in the crystallographic asymmetric unit. Each monomer has one THA molecule and two Mg²⁺ ions in the active center (Allen et al. 1995). We generated the Cartesian coordinates of the atoms of the second dimer because the tetramer is the active form of XyI. The structure of the natural substrate, D-xylose, was modeled in the four active centers by superimposing the common atoms between this molecule and THA. The O2 atom of D-xylose was modeled unprotonated. The short distance between the OD1 atom of Asp254 and the OE1 atom of Glu185 in the structure of the XyI–THA complex indicates a hydrogen bond interaction between them. Van Tilbeurgh et al. (1992) have proposed a mechanism for the O2 to O1 proton shuttle, in which the proton on O2 of D-xylose is transferred to Asp254 via the water ligand to Mg2. On the basis of these data, we decided that Asp254 was protonated in the active site. One of the active sites of the tetrameric enzyme was used to model the hydride transfer reaction, which was solvated with a 24 Å sphere of water molecules with the center coincide with the geometric center of atoms O2-C2-C1-O1 of the substrate. Water molecules that are less than 2.5 Å from any protein atom or crystallographic water were removed. The final system contained 25213 atoms. From this point on, the name 2GYI [PDB] will be used for this initial structure of the D-xylose-enzyme complex and for the structures obtained from MD simulations.

To model the Whitlow mechanism (Whitlow et al. 1991), we used the glucose–complex structure (PDB: 1XYB [PDB] ) as the starting coordinates to construct the reactive configuration for the hydride transfer reaction (Fig. 12, III). We first carried out PMF calculations with respect to the proton transfer reaction and its accompanying Mg–Mg migration (I III, Fig. 12). Then, the neutral Asp256 ligand of Mg2 was replaced by the carbonyl O1 atom of the substrate, and one of the bidentate oxygen ligands of Asp254 was displaced by the O2-alkoxide ion as a result of the proton transfer and Mg–Mg rearrangement. Note that these Mg–ligand exchanges involved no net loss or gain of charges, whereas a neutral water moleculate was proposed to replace the Asp254 anionic ligand (III, Fig. 12) in the mechanism of Allen et al. (1995).

In the XyI–glucose structure, the interaction between Asp254 and Glu185 is not present, and we rationalized that Asp256 was protonated, following the mechanism proposed by Whitlow et al. (1991). In the following discussion, this initial configuration and the structures obtained by running MD simulations started from it will be called 1XYB [PDB] .

The two models of the D-xylose-enzyme complex described above were heated and equilibrated by running stochastic boundary MD at 298 K (Garcia-Viloca et al. 2003), which provided the Michaelis complex structures used to start the calculations of the PMF by means of the umbrella sampling technique. More details of the model used have been explained in a recent article that contains the final results of our dynamical study (Garcia-Viloca et al. 2003). Here, we have the objective of elucidating structural issues by comparing the results obtained from the two different X-ray structures by using the same computational techniques and potential energy function applied to different initial enzyme-substrate configurations.

Next we turn attention to a comparison of the calculated PMFs and the interaction distances along the hydride transfer reaction. Figure 13 compares the PMFs obtained from the simulation based on the XyI–THA structure (2GYI [PDB] , solid line) and from the simulation based on the XyI–glucose structure (1XYB [PDB] , dotted line). The reaction coordinate, z, was defined as the difference between C2–H and C1–H distances (Fig. 8). Figure 13 indicates that the active site of the 2GYI [PDB] starting configuration provided weaker stabilization of the transition state and product complexes by 5 kcal/mole and 12 kcal/mole, respectively, relative to that in the second simulation.

View larger version (12K): [in this window] [in a new window]   Figure 13. Classical potentials of mean force obtained from the simulation based on the 2GYI [PDB] structure (solid line) and from the simulation based on the 1XYB [PDB] structure (dotted line). For details of these starting structures, see "XyI: The use of an inhibitor–enzyme complex proposed to be a transition state analog." By definition, each PMF is zero at the reactant state.

As stated above, both simulations were started from the proposed reactive conformation for the hydride transfer step (Collyer et al. 1990; Whitlow et al. 1991; Lavie et al. 1994; Allen et al. 1995; Garcia-Viloca et al. 2003), in which Mg2 has approached Mg1 after proton transfer from 2-hydroxyl group of the substrate to a hydroxide ion coordinated to Mg2. Accordingly, a similar Mg1–Mg2 average distance was obtained for the two simulations at the reactant state. However, the increase of this distance in going from the reactant to the product state is larger for the 1XYB [PDB] structures, and only in this case does it arrive to a value that corresponds to position 1 for Mg1. Thus, only the simulation based on the XyI–glucose structure (PDB 1XYB [PDB] ) is able to reproduce the full range of the breathing motion of the metal ions in the active site of XyI. The implications of this motion, which confirm the essential features of the mechanism postulated in whole or in part by Collyer et al. (1990) and Whitlow et al. (1991) and by Ringe, Petsko, and coworkers (Lavie et al. 1994; Allen et al. 1995), have been discussed in our previous paper (Garcia-Viloca et al. 2003).

A possible reason for the different behavior of Mg2 in the two structures may be inferred from the analysis of the interactions of Mg2 with its ligands. For the reactant, the distances between Mg2 and its ligands are very similar in the two sets of structures. However, an anionic ligand, Asp254, has been replaced by a neutral water molecule in the 2GYI [PDB] structures, which has a weaker interaction with the metal. This seems to have an effect on the distances between the metal and the O2 and O1 substrate atoms, which are shorter (1.97 Å and 1.99 Å, respectively) in the 2GYI [PDB] structures than in the 1XYB [PDB] structures (2.01 Å and 2.08 Å, respectively); this indicates a stronger interaction with the metal in the first case. The same trend is seen at the transition state. At the product state of the second simulation (1XYB [PDB] structures) Mg2 has lost its interaction with O2, which has been replaced by the interaction with a water molecule. However, the distance between the metal and O2 atom of D-xylose in the product structures obtained from the first simulation (2GYI [PDB] ) indicates that the interaction is still present, retaining Mg2 in a position closer to Mg1. The different coordination spheres of Mg2 in the two simulations are related to the different conformation of Asp254 in the active center of the initial X-ray structures and seem to be the main reason for the different results obtained.

In a previous article (Garcia-Viloca et al. 2003), we have presented the results of a model based on the 1XYB [PDB] structure. As mentioned in the Materials and Methods section, the energetic results obtained in that study are not identical to the ones obtained here for the model based on the same X-ray structure (PDB 1XYB [PDB] ), but the structural results obtained are comparable. These results suggest that the metal movement along the hydride transfer coordinate is reproduced by the simulations if Asp254 maintains one ligand to Mg2. In addition, the higher barrier obtained with the simulation based on the inhibitor structure (PDB 2GYI [PDB] ) does not support the loss of the two interactions between Mg2 and Asp254 after the proton transfer step prior to the hydride transfer step, which has been proposed by Petsko and coworkers (Lavie et al. 1994; Allen et al. 1995). Contrarily to the structural data found by these investigators in the structure of a transition state analog, in the two structures resolved by Collyer et al. (1990), which were also proposed to be analogs of the transition state, the OD2 atom of Asp254 is interacting with Mg2 (OD2–Mg2 distance is 2.4 Å and 2.7 Å).

Concluding remarks
In this article we have presented two case studies that show that small differences in the network of interactions around the active site can have profound influences on enzymatic reaction rates. In the first example, the study of the first step of the -oxidation cycle catalyzed by ACAD, two different crystal structures (obtained from different members of this class of enzymes) were used as starting points for simulations of the coupled proton transfer and hydride transfer. For both charge transfer steps, pronounced differences were found between the PMF calculated by using the two crystal structures. Specifically, for the proton transfer the free energy barrier differed by 13 kcal/mole, and for the hydride transfer the free energy barrier differed by 17 kcal/mole. The chief difference between the two structures is attributed to an interaction between the catalytic basic residue (Glu367 for BCAD and Glu376 for MCAD) and a nearby arginine (Arg247 for BCAD and Arg256 for MCAD). For one simulation, this interaction causes the species in the active center to be in a ion pair conformation that is not suitable for the proton transfer. Moreover, the substrate moves away from the FAD cofactor as a result of the proton abstraction, which also affects the simulation of the hydride transfer.

In the second example, the study of the hydride transfer reaction catalyzed by XyI, two Michaelis complex models were built by using two different X-ray structures of the enzyme from S. olivochromogenes. These two X-ray structures were resolved by the same investigators, under the same experimental conditions, but with a different ligand in the active center. Essentially, they only differ in the orientation of two residues in the active site. However, the results obtained by using the same methodology to generate configurations along the hydride transfer reaction path indicate that the relative energies obtained (from one or the other simulation) for the ensembles of reactant, transition state, and product structures are quite different. As a consequence, the two free energies of activation differ by 5 kcal/mole. The two simulations show a different trend in the changes of the interactions of one of the two Mg²⁺ cofactors with the other metal, with the substrate, and with the enzymatic residues during the hydride transfer. Our results indicate that the movement of this metallic cofactor during the reaction is very sensitive to the subtle differences in its coordination sphere in the two initial structures. This movement, which has been inferred from X-ray data, is found to be essential to lower the free energy of activation, and it is only reproduced in one of the two simulations.

These case studies illustrate the fundamental sensitivity of catalytic efficiency to even small differences in protein conformation near the active site, and they indicate that X-ray structures, because they are usually obtained with an empty active site or an inhibitor, may sometimes be unsuitable models for enzymatic simulations. One can take an even broader view and relate the present test cases to two long-standing and widely recognized questions: (1) To what extent does the structure in the crystal conform to the structure of the functional protein in a cell, and (2) to what extent can a MD simulation make up for a starting structure that happens to be in a nonproductive conformational local minimum? Although in principle a correct calculation of the PMF by umbrella sampling will be independent of the starting point, in practice, on practical time scales, the simulation will not explore regions separated by high barriers from the starting geometry unless those are explicitly removed by the umbrella potential. Furthermore, state-of-the art simulation techniques for enzyme reactions recognize that similar enzyme configurations can lead to very different reaction properties and that a quantitative rate calculation should include averaging over an ensemble of protein/substrate conformations (Alhambra et al. 2000; Gao and Truhlar 2002; Garcia-Viloca et al. 2002 Garcia-Viloca et al. 2003; Poulsen et al. 2003; Tresadern et al. 2003). Nevertheless, a critical time-scale issue remains. Functioning enzymes with reactive time scales of milliseconds explore all conformational states available through transitions that occur on the millisecond, microsecond, and faster time scales. But MD simulations are often carried out on time scales of a nanosecond or less; thus if the crystal structure used as a starting point for the simulation does not correspond to a typical productive structure of the enzyme, one may draw incorrect conclusions about the mechanism or underestimate the catalytic power of the enzyme.

Although the present article emphasized differences in X-ray structures, one would expect similar issues to arise if one starts out with structures obtained by molecular docking computations. Furthermore, when a mutant structure with substrate bound is used as the model for a Michaelis complex of the wild type by mutating the residues back, the potential changes in conformation caused by the substrate may be different in the wild type and the mutant. In fact, one might expect that a wide variety of starting structures might be available by combining available experimental data with a variety of computational protocols. As molecular modeling of protein function relies more and more heavily on computation, it is very important to be aware of the uncertainties in simulations of enzyme reactions that arise from the choice of starting structure. The many successes of MD simulations that have been reported along with the tendency (we imagine) to publish only one’s successful efforts have lessened the attention paid to this problem, but we hope that the two case studies presented here will be useful concrete examples of the seriousness of the problem.

	Footnotes

 
¹ Beinert 1963; Gopalan et al. 1984; Ikeda et al. 1985; Pohl et al. 1986; Finocchiaro et al. 1987; Powell and Thorpe et al. 1988; Schopfer et al. 1988; Bross et al. 1990; Kim et al. 1993; Ghisla et al. 1994; Shin et al. 1994; Djordjevic et al. 1995; Trivel et al. 1995; Dakoji et al. 1996; Lee et al. 1996; Nandy et al. 1996; Palfrey and Massey 1998; Rudik et al. 1998; Vock et al. 1998; Engst et al. 1999; Tamaoki et al. 1999; Pellett et al. 2000; Peterson et al. 2001; Bach et al. 2002; Gopalan and Srivastava 2002; Lamm et al. 2002; Poulsen et al. 2003; Ghisla and Thorpe 2004.

² Carrel et al. 1984; Farber et al. 1989; Collyer et al. 1990; Lee et al. 1990; Whitlow et al. 1991; Jenkins et al. 1992; Lambier et al. 1992; Rangarajan and Hartley 1992; van Tilbeurgh et al. 1992; Zheng et al. 1993; Allen et al. 1994a, b, 1995; Lavie et al. 1994; van Bastelaere et al. 1995; Bhosale et al. 1996; Nicoll et al. 2001; Garcia-Viloca et al. 2002, 2003.

	Acknowledgments

 
We are very grateful to Sandro Ghisla for sending his comments and those of J.J. Kim that helped us to improve an earlier version of this manuscript. This work was supported in part by a University of Minnesota Supercomputing Institute Research Scholarship to T.D.P. and by the NSF through a grant to D.G.T. T.D.P. also thanks Statens Naturvidenskabelige Forskningsråd for support.

The publication costs of this article were defrayed in part by payment of page charges. This article must therefore be hereby marked "advertisement" in accordance with 18 USC section 1734 solely to indicate this fact.

	References

TOP Abstract Introduction Materials and methods Results and Discussion References

 
Alhambra, C., Gao, J., Corchado, J.C., Villà, J., and Truhlar, D.G. 1999. Quantum mechanical dynamical effects in an enzyme-catalyzed proton-transfer reaction. J. Am. Chem. Soc. 121: 2253–2258.[CrossRef]

Alhambra, C., Corchado, J.C., Sanchez, M.L., Gao, J., and Truhlar, D.G. 2000. Quantum dynamics of hydride transfer in enzyme catalysis. J. Am. Chem. Soc. 122: 8197–8203.[CrossRef]

Allen, K.N., Lavie, A., Farber, G.K., Glasfeld, A., Petsko, G.A., and Ringe, D. 1994a. Isotopic exchange plus substrate and inhibition kinetics of D-xylose isomerase do not support a proton-transfer mechanism. Biochemistry 33: 1481–1487.