J. Phys. Chem. B 2010, 114, 3393–3397
3393
Importance of the Lactate Dehydrogenase Quaternary Structure in Theoretical Calculations Katarzyna S´widerek and Piotr Paneth* Institute of Applied Radiation Chemistry, Technical UniVersity of Lodz, ulica Zeromskiego 116, 90-924 Lodz, Poland ReceiVed: January 2, 2010; ReVised Manuscript ReceiVed: February 1, 2010
Using the example of lactate dehydrogenase, we show that enzyme quaternary structure has an important influence on the structure of the active site and that models that comprise all amino acids in the vicinity of an active site, but are missing this structural information, can lead to incorrect results. We also show that binding isotope effects are very sensitive to the geometric parameters, and thus one should be very cautious when interpreting results obtained with models that are too coarse. In terms of the type of hydrogen bonds, our results indicate that binding isotope effects are pronounced only when a hydrogen bond exhibits some covalent character. Introduction Enzymatic systems are very large, usually comprising many thousands of atoms, and therefore, they are not amenable to modeling at high levels of theory. The present remedy for this situation is a plethora of multilevel schemes, referred to in the literature as QM/MM methods.1 By embedding a quantum mechanics (QM) calculation in a classical molecular mechanics (MM) model of the environment, these methods are expected to capture the influence of the whole enzymatic environment on the reaction occurring in the enzyme active site. Although QM/MM methodology has rapidly developed in the past couple of years,2 the treatment of large systems is still limited by the computational cost and storage requirements. Recently, the size of the QM part;3 the performance of QM/QM methods in which the whole enzyme is treated quantum-mechanically;4 and approaches to calculations of very large, multimillion-atom systems have been discussed in literature.5 Because the active site of an enzyme is usually located between subunits, it is necessary to take into account both/all of these surrounding subunits, which significantly increases the size of the model. To facilitate calculations in such cases, distant parts of the enzyme are frequently skipped6 or frozen.7 Lactate dehydrogenase (LDH) provides an example of such a system: Based on MD simulations, it has been pointed out that a model that uses a single subunit is not sufficient and that two subunits should be used to describe the active site.8 We have carried out calculations of binding isotope effects (BIEs) for this enzyme and have shown that their values also depend on whether a monomer- or dimer-based model is used.9 In vitro LDH is a tetramer with two equivalent dimer units. One would therefore expect that the dimer should provide a sufficient model of the active site of this enzyme. The increase in computer power over recent years has allowed us to tackle this problem by comparing previous results obtained for monomer and dimer models with calculations carried out on the full tetrameric enzyme. As the starting point, we used the recently reported crystal structure [Protein Data Bank (PDB) code 3H3F]9 that contains both open (chain A) and closed (chain B) conformations of the active site. We have shown previously * Corresponding author. Phone: +48 42 631-3199. E-mail: paneth@ p.lodz.pl.
that B3LYP10/OPLS-AA11 implemented in the QSite program12 performed the best among the QM/MM and QM/QM approaches studied, and therefore, the same theoretical level was applied here to the full tetramer-based model. Computational Details The X-ray structure of LDH (PDB code 3H3F)9 was used as the starting point for modeling. All missing hydrogen atoms were added using the Maestro version 8.5 module of the Schrodinger software package.13 Additionally, His192 was protonated manually. Geometry optimization was subsequently performed using QM/MM schemes to model interactions of the inhibitor in the active site of models of LDH. First, the molecular mechanics force field OPLS-AA 200511 in combination with the B3LYP10 density functional theory functional expressed in the 6-31G++(d,p) basis set, as implemented in the QSite version 5.0 program,12 was used. To compare optimized structures, the standalone version of the Dali14 program (DaliLite15) was used with default parameters. The alignment of chain models obtained from the dimer and tetramer was achieved by rotation or translation of chain A or B to the position of the corresponding chain of monomers. These procedures were used interactively through the EMBLEBI (European Molecular Biology Laboratory European Bioinformatics Institute) Web page (http://www.ebi.ac.uk/Tools/ dalilite/index.html). The Ramachandran plots16 were obtained using the online Ramachandran plot 2.0 software hosted at the Supercomputer Education and Research Centre, Indian Institute of Science (http://dicsoft1.physics.iisc.ernet.in/rp/index.html). Theoretical values of binding isotope effects (BIEs) were subsequently calculated using the ISOEFF program.17 To obtain BIE values, oxamate in aqueous solution was modeled by placing an oxamate ion in a sphere of 516 water molecules. Using the HyperChem program,18 from a 20 × 20 × 20 Å cubic box of water molecules centered at the carboxylic carbon atom of oxamate, with a minimum distance between solvent and solute of 2.3 Å, a sphere with a radius of about 15 Å was prepared. Geometry optimization of this system was subsequently performed using the same QM/MM schemes as for the enzymatic system.
10.1021/jp100026z 2010 American Chemical Society Published on Web 02/15/2010
3394
S´widerek and Paneth
J. Phys. Chem. B, Vol. 114, No. 9, 2010
Figure 1. Superposed structures of chains A (left) and B (right) from monomer, dimer, and tetramer models.
Figure 2. Superposed active sites of monomer (yellow), dimer (blue), and tetramer (red) models.
Results and Discussion Structure Comparison. The most widely used statistic in protein structure comparisons is the root-mean-square distance (rmsd). rmsd measures the difference between the CR atom positions between structures of proteins. The smaller the deviation, the more spatially equivalent the structures. We define i)N δi2]1/2, where δ is rmsd for our studies as rmsd ) [(1/n)∑i)1 dimer/tetramer monomer - ri | between N pairs of the distance |ri equivalent CR atoms and ri represents the positions of the CR atoms. rmsd values were calculated using the DaliLite program. As the reference system, CR atom positions of the monomerbased models of chain A and B were used. The corresponding chains from the dimer and tetramer models were superimposed onto these models. Figure 1 presents two-dimensional plots illustrating the alignments of both types of chains with the calculated values of rmsd. The obtained results indicate that chains from theoretical models of different sizes, although not identical, are not substantially different. The largest deviation in the CR atom positions is observed in the part of the chains identified in Figure 1 as the “mobile tail”. In this contribution, we are focusing on changes that appear in the active sites of both chains (A and B) with open- and closed-loop conformations. For this reason, smaller models, containing oxamate and NADH molecules and the amino acid residues constituting the active site, were prepared. Their alignment is presented in Figure 2. It can be seen that there are meaningful differences in the active sites in chains A and B depending on the model. Thus, in chain A with the open-loop conformation, the largest differences between the positions of
the Gln99, Arg105, and Thr247 residues are observed, whereas in the active site of chain B, the largest changes concern the Arg105, His192, and Thr247 residues. To quantify these differences, the distances and torsion angles that change the most were studied in more detail. Selected distances to Gln99 and Arg105 of the open loop in chain A and hydrogen-bonding contacts with all amino acids constituting active sites obtained for different models are compared in Table 1, together with the crystallographic data. The atom numbering of the oxamate anion and its interactions with amino acids are given in Figure 3. As can be seen from the first three entries in Table 1, the distances from oxamate to the amino acids in the open loop (chain A: Glu99 and Arg105) are underestimated by all models, probably because of the high flexibility of the loop. The standard deviations of the calculated hydrogen-bonding contacts from the crystallographic results are equal to 0.5, 0.7, and 0.3 Å for the monomer, dimer, and tetramer, respectively. These results show that, even though the dimer-based model incorporates all amino acids surrounding the active sites, it does not adequately describe their geometries. Evidently, the quaternary structure plays an important role in tightening the active sites. Analysis of the hydrogen-bonding network in the active sites shows interesting variations between the models used. Table 2 lists distances, and Figure 4 presents the most prominent differences obtained for the studied models. Large values of H6-O3 and H3-O5 observed for Thr247 and Arg105 illustrate the deficiency of the monomer-based model. However, in chain B (the closed loop) of the dimer-based model, Glu99 is much farther from oxamate than in the tetramer-based model. The
Importance of LDH Quaternary Structure in Calculations
J. Phys. Chem. B, Vol. 114, No. 9, 2010 3395
TABLE 1: Comparison of Calculated Distances (Å) in the Active Sites with the Crystallographic Data chain A
chain B
distance
crystallography
monomer
dimer
tetramer
crystallography
monomer
dimer
tetramer
NGlu99-N6 NArg105-O5 NArg105-O1 NArg168-O1 NArg168-O3 NHis192-O5 OThr247-O3 CNADH-C1
16.23 13.18 15.81 2.79 3.45 3.76 2.53 3.12
10.11 8.82 9.77 2.82 2.95 3.39 3.8 3.7
6.16 9.74 10.16 2.9 2.84 3.46 2.79 4.08
8.21 8.98 9.66 2.84 2.85 3.29 2.73 3.7
3.4 2.96 3.08 2.87 3.25 3.07 3.29 3.35
3.64 3.22 2.85 2.83 2.86 3.21 3.48 3.62
5.21 3.16 2.89 2.83 2.93 2.92 2.89 4.78
3.51 3.07 2.92 2.89 2.96 3.06 2.96 3.73
same trend is observed for the distance between the C1 atom of oxamate and H7 of NADH, the hydrogen atom that is being transferred in the reaction with the physiological substrate, pyruvate. In the tetramer-based model, it is under 3 Å, whereas in the dimer-based model, it is nearly 5 Å! The above analysis indicates that the most important distance changes related to the changes in the hydrogen-bond strength concern the Thr247 residue. In the case of the monomer-based model, for both chains, the distance between this threonine and the oxamate molecule, after optimization, is longer than 4 Å. In the largest models, the H7-O3 distance decreases to about 2 Å. This distance shortening causes a significant increase of the hydrogenbond strength. The analysis of the spatial orientation of amino acid residues was carried out on the basis of a Ramachandran plot. Rotations in the case of an amino acid residue are possible around the N-CR (φ) and CR-C (ψ) single bonds. The Ramachandran plot shows the φ and ψ torsion angles. The sequence of the two mentioned angles and the angles of rotation about the peptide bond for all residues in the protein provide information regarding
Figure 3. Structure of active site of LDH with oxamate and NADH.
Figure 4. Major differences in distances between oxamate and activesite amino acids for different models.
the backbone conformation. The obtained results collected in Table 3 are best represented graphically, allowing for the identification of conformational differences of interesting amino acid residues between applied LDH. Ramachandran plots for the two chains, A and B, from the monomer, dimer, and tetramer models are presented in Figure 5. The values of the torsion angles confirm our observations that the largest deviations between the monomer-based and larger models are observed in chain A for Gln99, Arg105, and Thr247 and in chain B for Gln99 and Arg105. Deviations are not as large for chain B as they are for chain A. This can be explained by the closed conformation of the active site, which is tighter than in the case of the open-loop conformation, where the possibility of atom translocation is higher. Therefore, it is not surprising that significant differences in ψ and φ values for Gln99 and Arg105 in chain A were observed. However, in addition to residues that participate in the mobile loop, different torsion angles of the Thr247 residue also appear; in both chains, an increase of the φ angle and a TABLE 3: Torsion Angles ψ (°) and O (°) for Amino Acid Residues from the Active Site of LDH monomer
TABLE 2: Distances (Å) between Oxamate and Active-Site Amino Acids Obtained from Theoretical Calculations chain A
chain B
bond type monomer dimer tetramer monomer dimer tetramer H1-O5 H2-O1 H3-O5 H4-O1 H5-O3 H6-O3 H8-N6 H7-C1
8.03 8.94 2.82 1.80 1.93 4.65 9.14 3.13
9.17 9.75 2.97 1.87 1.81 1.83 5.93 3.74
8.23 8.89 2.75 1.81 1.82 1.77 7.93 3.00
2.33 1.84 3.61 1.82 1.83 4.06 3.64 2.89
2.20 1.86 1.93 1.80 1.90 1.96 5.05 4.82
2.12 1.92 2.15 1.87 1.93 2.02 3.41 2.91
residue
φ (°)
ψ (°)
Gln99 Arg105 Arg168 His192 Thr247
-14.6 -102.9 -70.9 -68.0 -116.1
-126.5 -175.5 -42.9 137.1 154.5
Gln99 Arg105 Arg168 His192 Thr247
-62.5 -52.4 -69.1 -68.5 -81.3
136.2 -37.8 -46.1 119.7 136.9
dimer φ (°)
ψ (°)
tetramer φ (°)
ψ (°)
Chain A -45.1 -68.0 -164.1 -163.3 -67.7 -42.8 -68.4 127.4 -131.3 142.3
-13.2 -162.5 -68.2 -78.0 -142.8
-114.7 -174.0 -45.0 138.1 113.0
Chain B -67.8 -63.8 -62.3 -66.9 -117.7
-73.1 -34.1 -67.0 -78.6 -118.8
135.8 -51.3 -48.1 120.5 127.1
135.0 -31.4 -52.4 128.5 127.4
3396
S´widerek and Paneth
J. Phys. Chem. B, Vol. 114, No. 9, 2010
Figure 5. Ramachandran plot of torsion angles ψ (°) and φ (°) for amino acid residues from the active site of LDH.
Figure 6. Binding isotope effects of oxamate heavy atoms: (left) Chain A (right)and chain B.
decrease of the ψ angle are observed. Deviations from the structure obtained for the monomer-based model in chain A are about 15° and 27° for φ and about 12° and 41° for ψ, for the dimer and tetramer models, respectively. In chain B, the differences between the monomer-based model and the dimer and tetramer models are the same and are equal to about 37° and 10° for angles φ and ψ, respectively. These changes in the torsion angle values, together with the previous analysis of distances in the active site, show that Thr247 is the residue that is most sensitive to changes of the model. Binding Isotope Effects. The most striking observation from the above discussion is the fact that, even though the dimerbased model includes all residues close to the active site, it differs from the larger model that is based on the full tetrameric form of the enzyme, and thus, it might not capture properly all important interactions within the active sites. We have addressed this problem by calculating BIEs, which are believed to be a sensitive probe of these interactions. Calculated BIEs of oxamate anion heavy atoms are presented in Figure 6. We start our discussion with the value of BIEs calculated for simultaneous isotopic substitution of both oxygen atoms (O1 and O3) of the carboxyl group (O1+O3BIE), as this is the only computed result that can be compared to the experimental value. We previously found19 its value for a singly labeled oxamate, which, on the basis of the rule of geometric mean,20 allows us to estimate the O1+O3BIE as about 0.97. As can be seen, the result obtained for chain A using the tetramerbased model is the closest to this value. Importantly, the values obtained for this BIE on the basis of the dimeric structure for the two chains differ only negligibly, whereas when the full tetrameric structure is included, these values differ by about 0.005, which makes them experimentally distinguishable by isotope-ratio mass spectrometry. On the contrary, nitrogen BIEs seem to be indicative of active-site conformation based on the dimer-based model; the expected value for chain A is positive (larger than unity) and inverse (less than unity) for chain B. However, the results obtained with the tetramer-based model show that this is not the case and that nitrogen BIEs are not useful in determining the active-site conformation.
Figure 7.
O3
BIE as a function of O3-H6 distance.
Binding isotope effects result from interactions of the reactants with the residues of the active site. The predominant interactions are hydrogen bonds. It is thus interesting to compare the magnitude of the observed BIEs with the strength of the hydrogen bonds. As indicated in Figure 3, O1+O3BIEs should reflect hydrogen bonds with Arg105, Arg168, and Thr247. As discussed above, of these three, the largest differences between the models are shown by Thr247. We have therefore studied how changes in the distance between the H6 proton from Thr247 and the O3 oxygen of oxamate (with the O3-H6-OThr247 angle fixed at 180°) affect the calculated O3BIEs. The results obtained show that the O3BIE starts to deviate from unity only when the O3-H6 distance becomes smaller than about 1.8 Å, at which point the O3BIE drops sharply, as illustrated in Figure 7. Interestingly, this threshold agrees very well with recent studies that indicated that this is the region where the covalent contribution starts to participate in a hydrogen bond21 and becomes prevalent below 1.2 Å. This supports the paradigm that isotope effects result from changes in force constants of reacting bonds and shows that hydrogen bonding is not an exception. Conclusions In summary, we have shown that the quaternary structure of LDH has a significant bearing on the structure of the active
Importance of LDH Quaternary Structure in Calculations site and that models that comprise all amino acids in the vicinity of the LDH active site but missing this structural information might lead to erroneous results. Because of limited resources and/or the lack of availability of a protein structure, in the modeling of enzymatic processes, putative residues of the active site4,22 or residues of the active site only6,23,24 are sometimes taken into account. Even when the whole enzyme embedded in water is considered, frequently the outermost sphere is frozen.7,25,26 Our findings indicate that these approaches can fail in some cases, such as the LDH binding studied herein. This observation, together with recently reported results suggesting that a few hundred atoms should be used in the QM part of QM/MM calculations,3 indicate that the modeling of enzymatic systems remains a challenging task. Second, we have shown that binding isotope effects are very sensitive to geometric parameters, and thus, one should be very cautious when interpreting results obtained with models that are too crude. Finally, in terms of the type of hydrogen bonds, our results indicate that binding isotope effects are pronounced only when a hydrogen bond exhibits some covalent character. Acknowledgment. This work was supported by Grant NN204/1579/33 from the Ministry of Science and Higher Education, Warsaw, Poland. Access to Polish national computation facilities at Cyfronet, Cracow; PSSC, Poznan; and ICM, Warsaw, is acknowledged. References and Notes (1) Warshel, A.; Levitt, M. J. Mol. Biol. 1976, 103, 227. (2) York, D. M., Lee, T.-S., Eds. Multi-scale Quantum Models for Biocatalysis: Modern Techniques and Applications; Challenges and Advances in Computational Chemistry and Physics Series; Leszczynski, J., Series Ed.; Springer: New York, 2009; Vol. 7. (3) (a) Lundberg, M; Siegbahn, P. E. M.; Morokuma, K. Biochemistry 2008, 47, 1031. (b) Solt, I.; Kulhanek, P.; Simon, I.; Winfield, S.; Payne, M. C.; Csanyi, G.; Fuxreiter, M. J. Phys. Chem. B 2009, 113, 5728. (4) Kwiecien, R.; Lewandowicz, A.; Paneth, P. Substrate-Enzyme Interactions from Modeling and Isotope Effects. In Molecular Materials with Specific Interactions: Modeling and Design; Sokalski, W. A., Ed.;
J. Phys. Chem. B, Vol. 114, No. 9, 2010 3397 Challenges and Advances in Computational Chemistry and Physics Series; Leszczynski, J., Series Ed.; Springer: New York, 2007; Vol. 4, pp 341363. (5) Schulz, R.; Lindner, B.; Petridis, L.; Smith, J. C. J. Chem. Theory Comput. 2009, 5, 2798. (6) (a) Mroginski, M. A.; Mark, F.; Thiel, W.; Hildebrandt, P. Biophys. J. 2007, 93, 1885. (b) Dyguda, E.; Grembecka, J.; Sokalski, W. A.; Leszczynski, J. J. Am. Chem. Soc. 2005, 127, 1658. (c) Dy´az, N.; Suarez, D.; Merz, K. M. J. Am. Chem. Soc. 2000, 122, 4197. (7) (a) Hart, J. C.; David, W.; Sheppard, D. W.; Hillier, I. H.; Burton, N. A. Chem. Commun. 1999, 79. (b) Ferrer, S.; Silla, E.; Tunon, I.; Oliva, M.; Moliner, V.; Williams, I. H. Chem. Commun. 2005, 5873. (c) Faulder, P. F.; Tresadern, G.; Chohan, K. K.; Scrutton, N. S.; Sutcliffe, M. J.; Hillier, I. H.; Burton, N. A. J. Am. Chem. Soc. 2001, 123, 8604. (8) Schmidt, R. K.; Gready, J. E. J. Mol. Model. 1999, 5, 153. (9) S´widerek, K.; Panczakiewicz, A.; Bujacz, A.; Bujacz, G.; Paneth, P. J. Phys. Chem. B 2009, 113, 12782. (10) (a) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (b) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (c) Miehlich, B.; Savin, A.; Stoll, H.; Preuss, H. Chem. Phys. Lett. 1989, 157, 200. (11) (a) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. 1996, 118, 11225. (b) Jorgensen, W. L.; McDonald, N. A. J. Mol. Struct. (THEOCHEM) 1998, 424, 145. (c) Jorgensen, W. L. BOSS, version 4.3; Yale University: New Haven, CT, 2001. (12) QSite, version 5.0; Schrodinger, LLC: New York, 2007. (13) Maestro, version 8.5; Schro¨dinger, LLC: New York, 2007. (14) Holm, L.; Park, J. Science 1996, 273, 595. (15) Holm, L.; Park, J. Bioinformatics 2000, 16, 566. (16) Ramanchadran, G. N.; Ramanchadran, C.; Sasisekharan, V. J. Mol. Biol. 1963, 7, 95. (17) Anisimov, V.; Paneth, P. J. Math. Chem. 1999, 26, 75. (18) HyperChem, release 8; Hypercube, Inc.: Gainesville, FL, 2007. (19) Gawlita, E.; Paneth, P.; Anderson, V. E. Biochemistry 1995, 34, 6050. (20) Bigeleisen, J. J. Chem. Phys. 1995, 23, 2264. (21) Grabowski, S. J.; Sokalski, W. A.; Dyguda, E.; Leszczynski, J. J. Phys. Chem. B 2006, 110, 6444. (22) Zhang, X.; DeChancie, J.; Gunaydin, H.; Chowdry, A. B.; Clemente, F. R.; Smith, A. J. T.; Handel, T. M.; Houk, K. N. J. Org. Chem. 2008, 73, 889. (23) Zhou, B.; Wong, C. F. J. Phys. Chem. A 2009, 113, 5144. (24) Idupulapati, N. B.; Mainardi, D. S. J. Phys. Chem A 2010, 114, 1887. (25) Dybala-Defratyka, A.; Paneth, P.; Banerjee, R.; Truhlar, D. G. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 10779. (26) Shaik, S.; Cohen, S.; Wang, Y.; Chen, H.; Kumar, D.; Thiel, W. Chem. ReV. 2010, 100, 949.
JP100026Z
