Thermodynamic Cycle Perturbation Research Articles

Catalytic Reaction Mechanism of Glyoxalase II: A Quantum Mechanics/Molecular Mechanics Study.

Methylglyoxal (MG) is a reactive and toxic compound produced in carbohydrate, lipid, and amino acid metabolism. The glyoxalase system is the main detoxifying route for MG and consists of two enzymes, glyoxalase I (GlxI) and glyoxalase II (GlxII). GlxI catalyzes the formation of S-d-lactoylglutathione from hemithioacetal, and GlxII converts this intermediate to d-lactate. A relationship between the glyoxalase system and some diseases like diabetes has been shown, and inhibiting enzymes of this system may be an effective means of controlling certain diseases. A detailed understanding of the reaction mechanism of an enzyme is essential to the rational design of competitive inhibitors. In this work, we use quantum mechanics/molecular mechanics (QM/MM) calculations and energy refinement utilizing the big-QM and QM/MM thermodynamic cycle perturbation methods to propose a mechanism for the GlxII reaction that starts with a nucleophilic attack of the bridging OH- group on the substrate. The coordination of the substrate to the Zn ions places its electrophilic center close to the hydroxide group, enabling the reaction to proceed. Our estimated reaction energies are in excellent agreement with experimental data, thus demonstrating the reliability of our approach and the proposed mechanism. Additionally, we examined alternative protonation states of Asp-29, Asp-58, Asp-134, and the bridging hydroxide ion in the catalytic process. However, these give less favorable reactions, a poorer reproduction of the crystal structure geometry of the active site, and higher root-mean-squared deviations of the active site residues in molecular dynamics simulations.

Reaction Mechanism of [NiFe] Hydrogenase Studied by Computational Methods.

[NiFe] hydrogenases catalyze the reversible conversion of molecular hydrogen to protons and electrons. This seemingly simple reaction has attracted much attention because of the prospective use of H2 as a clean fuel. In this paper, we have studied the full reaction mechanism of this enzyme with various computational methods. Geometries were obtained with combined quantum mechanical and molecular mechanics (QM/MM) calculations. To get more accurate energies and obtain a detailed account of the surroundings, we performed big-QM calculations with 819 atoms in the QM region. Moreover, QM/MM thermodynamic cycle perturbation calculations were performed to obtain free energies. Finally, density matrix renormalisation group complete active space self-consistent field calculations were carried out to study the electronic structures of the various states in the reaction mechanism. Our calculations indicate that the Ni-L state is not involved in the reaction mechanism. Instead, the Ni-C state is reduced by one electron and then the bridging hydride ion is transferred to the sulfur atom of Cys546 as a proton and the two electrons transfer to the Ni ion. This step turned out to be rate-determining with an energy barrier of 58 kJ/mol, which is consistent with the experimental rate of 750 ± 90 s-1 (corresponding to ∼52 kJ/mol). The cleavage of the H-H bond is facile with an energy barrier of 33 kJ/mol, according to our calculations. We also find that the reaction energies are sensitive to the size of the QM system, the basis set, and the density functional theory method, in agreement with previous studies.

Human Ferrochelatase: Insights for the Mechanism of Ferrous Iron Approaching Protoporphyrin IX by QM/MM and QTCP Free Energy Studies.

Ferrochelatase catalyzes the insertion of ferrous iron into protoporphyrin IX, the terminal step in heme biosynthesis. Some disputes in its mechanism remain unsolved, especially for human ferrochelatase. In this paper, high-level quantum mechanical/molecular mechanics (QM/MM) and free-energy studies were performed to address these controversial issues including the iron-binding site, the optimal reaction path, the substrate porphyrin distortion, and the presence of the sitting-atop (SAT) complex. Our results reveal that the ferrous iron is probably at the binding site coordinating with Met76, and His263 plays the role of proton acceptor. The rate-determining step is either the first proton removed by His263 or the proton transition within the porphyrin with an energy barrier of 14.99 or 14.87 kcal/mol by the quantum mechanical thermodynamic cycle perturbation (QTCP) calculations, respectively. The fast deprotonation step with the conservative residues rather than porphyrin deformation found in solution provides the driving force for biochelation. The SAT complex is not a necessity for the catalysis though it induces a modest distortion on the porphyrin ring.

Reduction Potentials and Acidity Constants of Mn Superoxide Dismutase Calculated by QM/MM Free‐Energy Methods

We used two theoretical methods to estimate reduction potentials and acidity constants in Mn superoxide dismutase (MnSOD), namely combined quantum mechanical and molecular mechanics (QM/MM) thermodynamic cycle perturbation (QTCP) and the QM/MM-PBSA approach. In the latter, QM/MM energies are combined with continuum solvation energies calculated by solving the Poisson-Boltzmann equation (PB) or by the generalised Born approach (GB) and non-polar solvation energies calculated from the solvent-exposed surface area. We show that using the QTCP method, we can obtain accurate and precise estimates of the proton-coupled reduction potential for MnSOD, 0.30±0.01 V, which compares favourably with experimental estimates of 0.26-0.40 V. However, the calculated potentials depend strongly on the DFT functional used: The B3LYP functional gives 0.6 V more positive potentials than the PBE functional. The QM/MM-PBSA approach leads to somewhat too high reduction potentials for the coupled reaction and the results depend on the solvation model used. For reactions involving a change in the net charge of the metal site, the corresponding results differ by up to 1.3 V or 24 pK(a) units, rendering the QM/MM-PBSA method useless to determine absolute potentials. However, it may still be useful to estimate relative shifts, although the QTCP method is expected to be more accurate.

An analysis of hydrophobic interactions of thymidylate synthase with methotrexate: Free energy calculations involving mutant and native structures bound to methotrexate

Since the human body for many reasons can adapt and become resistant to drugs, it is important to develop and validate computer aided drug design (CADD) methods that could help predict binding affinity changes that can result from these resistant enzymes. The free energy perturbation (FEP) methodology is the most accurate means of estimating relative binding affinities between inhibitors and protein variants. In this paper, we describe the role played by hydrophobic residues lining the active site region, particularly (79)Ile and (176)Phe, in the binding of methotrexate to the Escherichia coli (E. coli) thymidylate synthase (TS) enzyme, using the thermodynamic cycle perturbation (TCP) approach. The computed binding free energy differences on the binding of methotrexate to the native and some mutant E. coli TS structures have been compared with experimental results. Computationally, four different 'mutations' have been simulated on the TS enzyme with methotrexate (MTX): (79)Ile --> (79)Val; (79)Ile --> (79)Ala; (79)Ile --> (79)Leu; and (176)Phe --> (176)Ile. The calculated results indicate that in each of these cases, the native residues ((79) Ile and (176) Phe) interact more favorably with methotrexate than the mutant residues and these results are corroborated by experimental measurements. Binding preference to wild type residues can be rationalized in terms of their better hydrophobic contacts with the phenyl ring of methotrexate.

Accurate QM/MM Free Energy Calculations of Enzyme Reactions: Methylation by Catechol O-Methyltransferase.

We recently described a method to compute accurate quantum mechanical free energies [Rod, T. H.; Ryde, U. Phys. Rev. Lett. 2005, 94, 138302]. The method, which we term quantum mechanical thermodynamic cycle perturbation (QTCP), employs a molecular mechanics force field to sample phase space and, subsequently, a thermodynamic cycle to estimate QM/MM free energy changes. Here, we discuss the methodology in detail and test an approach based on a different thermodynamic cycle. We also show that a new way of treating hydrogen link atoms makes the free energy changes converge faster and that extrapolation to higher accuracy can be performed. We finally discuss the quantum mechanical free energy (QM/MM-FE) method in the framework of the QTCP method. All methods considered are applied to the methylation of catecholate catalyzed by catechol O-methyltransferase. We compute the free energy barrier for the reaction by computing free energy changes in steps between fixed QM regions along a predetermined reaction pathway. Using the QTCP approach, an extrapolated activation free energy of 69 kJ/mol for the forward reaction and 90 kJ/mol for the reverse reaction are obtained at the level of the B3LYP functional and the 6-311++G(2d,2p) basis set. The value for the forward reaction is in excellent agreement with the experimental value of 75 kJ/mol. Results based on the QM/MM-FE method differ by less than 10 kJ/mol from those values, indicating that QM/MM-FE may be a fairly accurate and cheap alternative to calculate QM/MM free energy changes. Moreover, the results are compared to barriers obtained with a fixed molecular mechanics environment as well as with structures optimized in a vacuum. All the computed free energy barriers are well converged. A major approximation in the current implementation of the QTCP method is that the QM region is fixed. The approximation leads to well-converged free energy barriers, which has been a problem in similar studies.

Solvation free energies calculated using the GB/SA model: Sensitivity of results on charge sets, protocols, and force fields

The sensitivity of aqueous solvation free energies (SFEs), estimated using the GB/SA continuum solvent model, on charge sets, protocols, and force fields, was studied. Simple energy calculations using the GB/SA solvent model were performed on 11 monofunctional organic compounds. Results indicate that calculated SFEs are strongly dependent on the charge sets. Charges derived from electrostatic potential fitting to high level ab initio wave functions using the CHELPG procedure and “class IV” charges from AM1/CM1a or PM3/CM1p calculations yielded better results than the corresponding Mulliken charges. Calculated SFEs were similar to MC/FEP energies obtained in the presence of explicit TIP4P water. Further improvements were obtained by using GVB/6-31G** and MP2/6-31+G** (CHELPG) charge sets that included correlation effects. SFEs calculated using charge sets assigned by the OPLSA* force field gave the best results of all standard force fields (MM2*, MM3*, MMFF, AMBER*, and OPLSA*) implemented in MacroModel. Comparison of relative and absolute SFEs computed using either the GB/SA continuum model or MC/FEP calculations in the presence of explicit TIP4P water showed that, in general, relative SFEs can be estimated with greater accuracy. A second set of 20 mono- and difunctional molecules was also studied and relative SFEs estimated using energy minimization and thermodynamic cycle perturbation (TCP) protocols. SFEs calculated from TCP calculations using the GB/SA model were sensitive to bond lengths of dummy bonds (i.e., bonds involving dummy atoms). In such cases, keeping the bond lengths of dummy bonds close to the corresponding bond lengths of the starting structures improved the agreement of TCP-calculated SFEs with energy minimization results. Overall, these results indicate that GB/SA solvation free energy estimates from simple energy minimization calculations are of similar accuracy and value to those obtained using more elaborate TCP protocols. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 769–780, 1998

An Approach to Rapid Estimation of Relative Binding Affinities of Enzyme Inhibitors: Application to Peptidomimetic Inhibitors of the Human Immunodeficiency Virus Type 1 Protease

This report describes a method for rapid assessment of the binding affinities of a series of analogous ligands to an enzyme. This approach is based on two variables (scores), representing (i) the enthalpy of binding and (ii) the strength of hydrophobic interaction. The method is then used to evaluate the binding of 11 different peptidomimetic inhibitors to the HIV-1 protease. Three-dimensional structures of these enzyme-inhibitor complexes are modeled based on the crystal structures of HIV-1 protease complexes with the known inhibitors. These structures are minimized using the AMBER force field, and the scores of binding enthalpy for each of the ligands are calculated. A second score to represent the hydrophobic interaction between a pair of atoms uses an exponential function of distance between the atoms and the product of their atomic hydrophobicity constants. This exponential function is used to assess the hydrophobic interaction energy between an enzyme and its inhibitor and also to compute and display a 'molecular hydrophobicity map' as a 3D visualization tool. These methods are then applied to obtain trends in relative binding affinities of pairs of analogous inhibitors. Calculated scores agree well with corresponding results from thermodynamic cycle perturbation (TCP) simulations as well as experimental binding data. Since the proposed calculations are computationally cheaper and faster than TCP calculations, it is suggested that these scores can form the basis for rapid, preliminary theoretical screening of proposed derivatives of an inhibitor prior to TCP analysis, synthesis, and testing.

Calculation of relative free energy differences for the covalent hydration of organic compounds: A combined quantum mechanical and free energy perturbation study

AbstractThe relative free energy difference (ΔΔGhyd) for the reversible addition of water to two unsaturated molecules is accurately computed using a combination of ab initio quantum mechanical calculations and free energy perturbation methods. Initial attempts to calculate the absolute hydration free energy difference (ΔGhyd) for formaldehyde and trichloroacetaldehyde gave values that differed substantially from experimental results even after inclusion of electron correlation energy contributions using third‐order (MP3) and fourth‐order (MP4) Møller‐Plesset perturbation theory and QCISD(T) correlation methods at the 6‐31G** basis set level. Inaccuracies in ΔGhyd were attributed to errors in the calculation of both ΔGgas and ΔΔGsol. Gas phase quantum mechanical free energies (ΔGgas) varied significantly (2–3 kcal/mol) depending on the level of theory. Errors in ΔΔGsol were attributed to slow convergence of the calculations using the thermodynamic cycle perturbation (TCP) method with explicit solvent. These errors were minimized or canceled, however, when relative hydration free energy differences (ΔΔGhyd) were calculated using a combination of ab initio quantum mechanical calculations and free energy perturbation methods. Calculated values for a variety of aldehydes and ketones were consistent with experimental data. © 1995 John Wiley & Sons, Inc.

Crystal-structure-based design and synthesis of novel C-terminal inhibitors of HIV protease.

The X-ray crystal-structure-based design, synthesis, computational evaluation, and activity of a novel class of HIV protease inhibitors are described. The initial lead compounds 2 and 3 were designed by modeling replacement groups for the C-terminal Val-Val-OCH3 of a known hydroxyethylene inhibitor into the active site of the reported crystal structure of HIV protease complexed with MVT-101. The lead compound 2 was found to be a modest inhibitor with a Ki = 1.67 microM. The X-ray crystal structure of compound 2 complexed with HIV protease was solved and used for subsequent design. The lead compound 3 was found to be a more potent inhibitor with Ki = 0.2 microM, and the structure of it complexed with HIV protease was also solved and used for subsequent design. Modification of both the C-terminus and N-terminus of indole 3 resulted in compounds with Ki = 30 nM. Using the crystal structure of compounds 2 and 3 with HIV protease as a starting point, the thermodynamic cycle perturbation molecular dynamics method was applied to a select group of compounds in order to test the accuracy of this type of computation within a series of closely related compounds.

Molecular recognition: effect of rotational isomers on host-guest binding

Monte Carlo simulations have been used to study the relative binding of dimethylurea and imidazolidone to a synthetic host molecule in CHCl 3 . The thermodynamic cycle-perturbation method was used to calculate the realtive free energy of binding, which was compared with experimental data from NMR binding studies. Special techniques have been used to properly account for the different rotational isomeric states of dimethylurea in the thermodynamic averages. The computed relative free energy of binding ΔΔG=3.6 kcal mol -1 facors the binding of imidazoline and compares reasonably well with the experimental value of 3.1 kcal mol -1

Calculation of solvation and binding free energy differences for folate-based inhibitors of the enzyme thymidylate synthase

The thermodynamic cycle perturbation approach, in conjunction with molecular dynamics simulations, has been used to calculate relative binding free energies for the closely related inhibitors 10-propargyl-5,8-dideazafolic acid (PDDF) and 10-formyl-5,8-dideazafolic acid (FDDF) to the binary complex consisting of the enzyme E. coli thymidylate synthase and 5-fluoro-2'-deoxyuridylate (FdUMP). The calculated difference in binding free energy (2.9 kcal/mol) is in good agreement with the experimentally observed preference for PDDF (3.75 kcal/mol)

Relative differences in the binding free energies of human immunodeficiency virus 1 protease inhibitors: a thermodynamic cycle-perturbation approach.

Peptidomimetic inhibitors of the human immunodeficiency virus 1 protease show considerable promise for treatment of AIDS. We have, therefore, been seeking computer-assisted drug design methods to aid in the systematic design of such inhibitors from a lead compound. Here we report thermodynamic cycle-perturbation calculations (using molecular dynamics simulations) to compute the relative difference in free energy of binding that results when one entire residue (valine) is deleted from one such inhibitor. In particular, we studied the "alchemic" mutation of the inhibitor Ac-Ser-Leu-Asn-(Phe-Hea-Pro)-Ile-Val-OMe (S1) to Ac-Ser-Leu-Asn-(Phe-Hea-Pro)-Ile-OMe (S2), where Hea is hydroxyethylamine, in two different (R and S) diastereomeric configurations of the hydroxyethylene group. The calculated (averaged for R and S) difference in binding free energy [3.3 +/- 1.1 kcal/mol (mean +/- SD); 1 cal = 4.184 J] is in good agreement with the experimental value of 3.8 +/- 1.3 kcal/mol, obtained from the measured Ki values for an equilibrium mixture of R and S configurations. Precise testing of our predictions will be possible when binding data become available for the two disastereomers separately. The observed binding preference for S1 is explained by the stronger ligand-protein interaction, which dominates an opposing contribution arising from the large desolvation penalty of S1 relative to S2. This calculation suggests that the thermodynamic cycle-perturbation approach can be useful even when a relatively large change in the ligand is simulated and supports the use of the thermodynamic cycle-perturbation algorithm for screening proposed derivatives of a lead inhibitor/drug prior to their synthesis.

Free Energy Calculations Predict Sequence Specificity In DNA-drug Complexes

Abstract Molecular dynamics simulations of DNA-netropsin complexes in water were performed using the thermodynamic cycle-perturbation method to calculate the free energy difference between complexes with an adenine-containing binding site and corresponding complexes where adenines are replaced by 2,6-diaminopurines (dap). The calculations predict a free energy difference of 3.7±0.9 kcal/mol (at 300K) in favour of netropsin binding to an (AATT)2 DNA sequence compared to a (dapdapTT)2, sequence.

Application of the thermodynamic cycle perturbation method to systems with nonadditive potentials

We have extended the thermodynamic perturbation cycle approach to compute thermodynamic properties of systems with nonadditive potentials. A rapid iterative method was developed and tested by application to systems of water microclusters of alkali metal cations and halogen anions interacting by means of a model potential incorporating polarization. The perturbation cycle calculations are in agreement with the results of direct calculation as long as both initial and final microclusters are unambiguously stable. For marginally stable cationic microclusters, e.g. Rb(H2O)+6 or Cs(H2O)+5, the perturbation cycle approach appears to provide more reliable estimates of thermodynamic properties. Transmutation entropy estimates entail substantial uncertainty, even for transitions between stable microclusters. They are in fair agreement with those of other theoretical studies but differ noticeably from the experimental values.

Thermodynamic cycle‐perturbation study of the binding of trifluoroacetyl dipeptide anilide inhibitors with porcine pancreatic elastase

The variety of results of crystallographic studies of the serine proteases complexed with isocoumarin inhibitors presents a challenging problem to modeling methods and molecular energetics. Therefore, the thermodynamic cycle-perturbation technique has been used to study a model system of elastase and two peptidic inhibitors. Using the program AMBER, the technique correctly predicts changes of the binding constants for the trifluoroacetyl dipeptide inhibitors in comparison with available experimental (kinetic and crystallographic) data. However, the absolute values obtained are shown to be sensitive to the specific electrostatic interaction potential parameters used in the simulations. The reader and user are cautioned that thermodynamic cycle-perturbation results may be too optimistic by underestimating the accuracy of free energy values. This is especially a matter of concern for those cases where a direct comparison with experimental values is not possible, viz., (1) the stimulation of binding of novel compounds, (2) structurally uncertain binding sites, or (3) structurally different binding modes. With our best 4-31G* ESP (electrostatic potential) charges we were able to reproduce experimentally determined free energy differences (delta delta A) with an accuracy of about 1.5 kcal/mol. Dynamically induced structural changes in the binding site of elastase, and particularly changes in hydrogen-bond patterns of the binding site, are also reported.

Computer-aided molecular design.

Theoretical chemistry, as implemented on fast computers, is beginning to yield accurate predictions of the thermodynamic and kinetic properties of large molecular assemblies. In addition to providing detailed insights into the origins of molecular activity, theoretical calculations can be used to design new molecules with specific properties. This article describes two types of calculations that show special promise as design tools, the thermodynamic cycle-perturbation method and the Brownian reactive dynamics method. These methods can be applied to calculate equilibrium and rate constants that describe many aspects of molecular recognition, stability, and reactivity.

Geometric considerations in the calculation of relative free energies of activation

A method that locates transition state structures between homologous reactions and, with the use of the thermodynamic cycle-perturbation technique, determines the relative free energy of activation between the transition states is presented. A simple model system which displays the problem of finding condensed phase transition states is used to illustrate the method.

Computer Method Improves Estimates of Molecular Properties

Using something of a thermodynamic sleight of hand, chemists at the University of Houston have developed a new method to estimate certain key properties of new molecules before the molecules are synthesized. J. Andrew McCammon, who is M. D. Anderson Professor of Chemistry at the university, and postdoctoral fellows Terry P. Lybrand and Chung F. Wong call their technique the thermodynamic cycle-perturbation approach. So far, they have used it to estimate the binding affinities for several types of receptor- binding interactions. In theory, though, the technique probably is much more widely applicable, perhaps to any chemical property that is based on thermodynamics. The need for a theoretical method that can accurately predict properties like binding affinity is great, McCammon points out. In pharmacology, for example, efforts to improve existing drugs or drug candidates often boil down to trying to find a modification that will improve the drug's binding to its target receptor or diminish it ...