Combinedab initioand Free Energy Calculations To Study Reactions in Enzymes and Solution: Amide Hydrolysis in Trypsin and Aqueous Solution

Integration of metal organic frameworks with enzymes as multifunctional solids for cascade catalysis

Enzymes exhibit a large degree of compatibility with metal-organic frameworks (MOFs) which allows the development of multicomponent catalysts consisting of enzymes adsorbed or occluded by MOFs. The combination of enzymes and MOFs in a multicomponent catalyst can be used to promote cascade reactions in which two or more individual reactions are performed in a single step. Cascade reactions take place due to the cooperation of active sites present on the MOF with the enzyme. A survey of the available data establishes that often an enzyme undergoes stabilization by association with a MOF and the system exhibits notable recyclability. In addition, the existence of synergism is observed as a consequence of the close proximity of all the required active sites in the multicomponent catalyst. After an introductory section describing the specific features and properties of enzyme-MOF assemblies, the main part of the present review focuses on the description of the cascade reactions that have been reported with commercial enzymes associated with MOFs, paying special attention to the advantages derived from the multicomponent catalyst. Related to the catalytic activity to metabolize glucose, generating reactive oxygen species (ROS) and decreasing the solution pH, an independent section describes the recent use of enzyme-MOF catalysts in cancer therapy. The last paragraphs summarize the current state of the art and provide our view on future developments in this field.

Efficient Computation of Free Energy Surfaces of Diels⁻Alder Reactions in Explicit Solvent at Ab Initio QM/MM Level

For Diels–Alder (DA) reactions in solution, an accurate and converged free energy (FE) surface at ab initio (ai) quantum mechanical/molecular mechanical (QM/MM) level is imperative for the understanding of reaction mechanism. However, this computation is still far too expensive. In a previous work, we proposed a new method termed MBAR+wTP, with which the computation of the ai FE profile can be accelerated by several orders of magnitude via a three-step procedure: (I) an umbrella sampling (US) using a semi-empirical (SE) QM/MM Hamiltonian is performed; (II) the FE profile is generated using the Multistate Bennett Acceptance Ratio (MBAR) analysis; and (III) a weighted Thermodynamic Perturbation (wTP) from the SE Hamiltonian to the ai Hamiltonian is performed to obtain the ai QM/MM FE profile using weight factors from the MBAR analysis. In this work, this method is extended to the calculations of two-dimensional FE surfaces of two Diels–Alder reactions of cyclopentadiene with either acrylonitrile or 1-4-naphthoquinone at ai QM/MM level. The accurate activation free energies at the ai QM/MM level, which are much closer to the experimental measurements than those calculated by other methods, indicate that this MBAR+wTP method can be applied in the studies of complex reactions in condensed phase with much-enhanced efficiency.

QM/MM methods for free energies and photochemistry

Hybrid computational methods describing a small region of a biomolecular system with quantum mechanics and the bulk with molecular mechanics, referred to as QM/MM methods, are now a central part of computational biochemistry. This review considers developments in the QM/MM approach that make it easier to calculate free energies using accurate QM-based potential energy expressions. We also describe techniques to treat electronic coupling between the core region and the MM environment. Polarizability of the protein matrix is important but so is electronic coupling. Applications of these new methods, especially to photochemistry, are discussed.

Averaged Condensed Phase Model for Simulating Molecules in Complex Environments

The need for configurational sampling dramatically increases the cost of combined quantum mechanics/molecular mechanics (QM/MM) simulations of chemical processes in solution. We developed an averaged condensed phase environment (ACPE) model that constructs an effective polarizable environment directly from explicitly sampled molecular dynamics configurations via the K-means++ algorithm and a mathematically rigorous translation of the molecular mechanics parameters. The model captures detailed heterogeneous features in the environment that may be difficult to describe using a conventional polarizable continuum model. Instead of performing repeated QM/MM calculations for each new configuration of the environment, the ACPE approach allows one to perform a single QM calculation on an averaged configuration. Here, we demonstrate the model by computing electronic excitation energies for several small molecules in solution. The ACPE model predicts the excitation energies in excellent agreement with conventional configurational averaging yet with orders of magnitude of reduction in the computational cost.

Methyltransferases do not work by compression, cratic, or desolvation effects, but by electrostatic preorganization

The enzyme catechol O-methyltransferase (COMT) catalyzes the transfer of a methyl group from S-adenosylmethionine to dopamine and related catechols. The search for the origin of COMT catalysis has led to different proposals and hypothesis, including the entropic, the NAC, and the compression proposals as well as the more reasonable electrostatic idea. Thus, it is important to understand the catalytic power of this enzyme and to examine the validity of different proposals and in particular the repeated recent implication of the compression idea. The corresponding analysis should be done by well-defined physically-based considerations that involve computations rather than circular interpretations of experimental results. Thus, we explore here the origin of the catalytic efficiency of COMT by using the empirical valence bond and the linear response approximation approaches. The results demonstrate that the catalytic effect of COMT is mainly due to electrostatic preorganization effects. It is also shown that the compression, NAC and entropic proposals do not account for the catalytic effect.

Quantum mechanics/molecular mechanics dual Hamiltonian free energy perturbation

The dual Hamiltonian free energy perturbation (DH-FEP) method is designed for accurate and efficient evaluation of the free energy profile of chemical reactions in quantum mechanical/molecular mechanical (QM/MM) calculations. In contrast to existing QM/MM FEP variants, the QM region is not kept frozen during sampling, but all degrees of freedom except for the reaction coordinate are sampled. In the DH-FEP scheme, the sampling is done by semiempirical QM/MM molecular dynamics (MD), while the perturbation energy differences are evaluated from high-level QM/MM single-point calculations at regular intervals, skipping a pre-defined number of MD sampling steps. After validating our method using an analytic model potential with an exactly known solution, we report a QM/MM DH-FEP study of the enzymatic reaction catalyzed by chorismate mutase. We suggest guidelines for QM/MM DH-FEP calculations and default values for the required computational parameters. In the case of chorismate mutase, we apply the DH-FEP approach in combination with a single one-dimensional reaction coordinate and with a two-dimensional collective coordinate (two individual distances), with superior results for the latter choice.

Presence or absence of a novel charge-transfer complex in the base-catalyzed hydrolysis of N-ethylbenzamide or ethyl benzoate

Reaction paths of base-catalyzed hydrolyses of isoelectronic substrates, Ph–C(=O)–X–Et [X = O (ethyl benzoate) and X = NH (N-ethylbenzamide)], were traced by DFT calculations. To simulate bond interchanges accompanied by proton transfers, a cluster model of Ph–C(=O)–X–Et + OH⁻(H₂O)₁₆ was employed. For X = O, three elementary processes and for X = NH four ones were obtained. The rate-determining step of X = O is the first TS (TS1, the OH⁻ addition step), while that of X = NH is TS2. TS2 of X = NH leads to a novel Mulliken charge-transfer complex, Ph–(OH)(O=)C∙∙∙N(H₂)–Et. The superiority or inferiority between the direct nucleophilic process or the general base-catalyzed process for TS1 was examined with the model Ph–C(=O)–X–Et + OH⁻(H₂O)_n, n = 3, 5, 8, 12, 16, 24 and 32. The latter process was calculated to be more favorable regardless of the number (n, except n = 3) of water molecules. The counter ion Na⁺ works unfavorably on the ester hydrolysis, particularly on TS1. A minimal model of TS1 was proposed and was found to be insensitive to n.

A molecular mechanism of hydrolysis of peptide bonds at neutral pH using a model compound

The covalent stability of peptide bonds is a critical aspect of biological chemistry and therapeutic protein applications. In this computational study, the hydrolytic reaction of peptide bonds at neutral pH was studied using a model compound, N-MAA. The most probable reaction pathway and intermediate(s) involved are controversial in previous studies. In addition, most previous computational studies focus on the energetics of chemical species involved, rather than providing a dynamic picture of the reaction process in aqueous conditions. However, fluctuations at finite temperatures are quite important, as we show. Thus, a path sampling method was used to generate an ensemble of trajectories according to their statistical weights in trajectory space. An ab initio molecular dynamics technique was applied to advance the time of the reaction in order to collect trajectories. The likelihood maximization procedure and its modification were used in extracting dynamically relevant degrees of freedom in the system, and approximations of the reaction coordinate were compared. It was found that this hydrolytic reaction is very complex because it involves many degrees of freedom. The reaction coordinate C-O distance previously assumed was found to be inadequate in describing the dynamic progress of the reaction. In addition to affecting atoms directly involved in bond-making and -breaking processes, the water network also has determining effects on the hydrolytic reaction, a fact which is manifest in the expression of the best one-dimensional reaction coordinate that we found, which includes five geometric quantities. p(B) histograms were computed to verify the results of the likelihood maximization and to evaluate the accuracy of our best reaction coordinate to the "true" reaction coordinate. The relation with previous suggested reaction pathways and intermediate(s) is discussed in terms of computational system, method, and accuracy.

Absolute binding free energy calculations: on the accuracy of computational scoring of protein-ligand interactions

Calculating the absolute binding free energies is a challenging task. Reliable estimates of binding free energies should provide a guide for rational drug design. It should also provide us with deeper understanding of the correlation between protein structure and its function. Further applications may include identifying novel molecular scaffolds and optimizing lead compounds in computer-aided drug design. Available options to evaluate the absolute binding free energies range from the rigorous but expensive free energy perturbation to the microscopic linear response approximation (LRA/beta version) and related approaches including the linear interaction energy (LIE) to the more approximated and considerably faster scaled protein dipoles Langevin dipoles (PDLD/S-LRA version) as well as the less rigorous molecular mechanics Poisson-Boltzmann/surface area (MM/PBSA) and generalized born/surface area (MM/GBSA) to the less accurate scoring functions. There is a need for an assessment of the performance of different approaches in terms of computer time and reliability. We present a comparative study of the LRA/beta, the LIE, the PDLD/S-LRA/beta, and the more widely used MM/PBSA and assess their abilities to estimate the absolute binding energies. The LRA and LIE methods perform reasonably well but require specialized parameterization for the nonelectrostatic term. The PDLD/S-LRA/beta performs effectively without the need of reparameterization. Our assessment of the MM/PBSA is less optimistic. This approach appears to provide erroneous estimates of the absolute binding energies because of its incorrect entropies and the problematic treatment of electrostatic energies. Overall, the PDLD/S-LRA/beta appears to offer an appealing option for the final stages of massive screening approaches.

Molecular mechanism of acid-catalyzed hydrolysis of peptide bonds using a model compound

The stability of peptide bonds is a critical aspect of biological chemistry and therapeutic protein applications. Recent studies found elevated nonenzymatic hydrolysis in the hinge region of antibody molecules, but no mechanism was identified. As a first step in providing a mechanistic interpretation, this computational study examines the rate-determining step of the hydrolytic reaction of a peptide bond under acidic pH by a path sampling technique using a model compound N-MAA. Most previous computational studies did not include explicit water molecules, whose effects are significant in solution chemistry, nor did they provide a dynamic picture for the reaction process in aqueous conditions. Because no single trajectory can be used to describe the reaction dynamics due to fluctuations at finite temperatures, a variant version of the transition path sampling technique, the aimless shooting algorithm, was used to sample dynamic trajectories and to generate an ensemble of transition trajectories according to their statistical weights in the trajectory space. Each trajectory was computed as the time evolution of the molecular system using the Car-Parrinello molecular dynamics technique. The likelihood maximization procedure and its modification were used in extracting dynamically relevant degrees of freedom in the system, and approximations of the reaction coordinate were compared. Its low log-likelihood score and poor p(B) histogram indicate that the C-O distance previously assumed as the reaction coordinate for the rate-determining step is inadequate in describing the dynamics of the reaction. More than one order parameter in a candidate set including millions of geometric quantities was required to produce a convergent reaction coordinate model; its involvement of many degrees of freedom suggests that this hydrolytic reaction step is very complex. In addition to affecting atoms directly involved in bond-making and -breaking processes, the water network also has determining effects on the hydrolytic reaction, a fact that is manifest in the expression of the one-dimensional best-ranked reaction coordinate, which includes three geometric quantities. The p(B) histograms were computed to verify the results of the likelihood maximization and to verify the accuracy of approximation to the "true" reaction coordinate.

Combined QM/MM mechanistic study of the acylation process in furin complexed with the H5N1 avian influenza virus hemagglutinin's cleavage site

Combined quantum mechanical/molecular mechanical (QM/MM) techniques have been applied to investigate the detailed reaction mechanism of the first step of the acylation process by furin in which the cleavage site of the highly pathogenic avian influenza virus subtype H5N1 (HPH5) acts as its substrate. The energy profile shows a simultaneous mechanism, known as a concerted reaction, of the two subprocesses: the proton transfer from Ser368 to His194 and the nucleophilic attack on the carbonyl carbon of the scissile peptide of the HPH5 cleavage site with a formation of tetrahedral intermediate (INT). The calculated energy barrier for this reaction is 16.2 kcal.mol(-1) at QM/MM B3LYP/6-31+G*//PM3-CHARMM22 level of theory. Once the reaction proceeds, the ordering of the electrostatic stabilization by protein environment is of the enzyme-substrate < transition state < INT complexes. Asp153 was found to play the most important role in the enzymatic reaction by providing the highest degree of intermediate complex stabilization. In addition, the negatively charged carbonyl oxygen of INT is well stabilized by the oxyanion hole constructed by Asn295's carboxamide and Ser368's backbone.

QM/MM methods for biomolecular systems

Combined quantum-mechanics/molecular-mechanics (QM/MM) approaches have become the method of choice for modeling reactions in biomolecular systems. Quantum-mechanical (QM) methods are required for describing chemical reactions and other electronic processes, such as charge transfer or electronic excitation. However, QM methods are restricted to systems of up to a few hundred atoms. However, the size and conformational complexity of biopolymers calls for methods capable of treating up to several 100,000 atoms and allowing for simulations over time scales of tens of nanoseconds. This is achieved by highly efficient, force-field-based molecular mechanics (MM) methods. Thus to model large biomolecules the logical approach is to combine the two techniques and to use a QM method for the chemically active region (e.g., substrates and co-factors in an enzymatic reaction) and an MM treatment for the surroundings (e.g., protein and solvent). The resulting schemes are commonly referred to as combined or hybrid QM/MM methods. They enable the modeling of reactive biomolecular systems at a reasonable computational effort while providing the necessary accuracy.

Progress in ab initio QM/MM free-energy simulations of electrostatic energies in proteins: accelerated QM/MM studies of pKa, redox reactions and solvation free energies

Hybrid quantum mechanical/molecular mechanical (QM/MM) approaches have been used to provide a general scheme for chemical reactions in proteins. However, such approaches still present a major challenge to computational chemists, not only because of the need for very large computer time in order to evaluate the QM energy but also because of the need for proper computational sampling. This review focuses on the sampling issue in QM/MM evaluations of electrostatic energies in proteins. We chose this example since electrostatic energies play a major role in controlling the function of proteins and are key to the structure-function correlation of biological molecules. Thus, the correct treatment of electrostatics is essential for the accurate simulation of biological systems. Although we will be presenting different types of QM/MM calculations of electrostatic energies (and related properties) here, our focus will be on pKa calculations. This reflects the fact that pKa's of ionizable groups in proteins provide one of the most direct benchmarks for the accuracy of electrostatic models of macromolecules. While pKa calculations by semimacroscopic models have given reasonable results in many cases, existing attempts to perform pKa calculations using QM/MM-FEP have led to discrepancies between calculated and experimental values. In this work, we accelerate our QM/MM calculations using an updated mean charge distribution and a classical reference potential. We examine both a surface residue (Asp3) of the bovine pancreatic trypsin inhibitor and a residue buried in a hydrophobic pocket (Lys102) of the T4-lysozyme mutant. We demonstrate that, by using this approach, we are able to reproduce the relevant side chain pKa's with an accuracy of 3 kcal/mol. This is well within the 7 kcal/mol energy difference observed in studies of enzymatic catalysis, and is thus sufficient accuracy to determine the main contributions to the catalytic energies of enzymes. We also provide an overall perspective of the potential of QM/MM calculations in general evaluations of electrostatic free energies, pointing out that our approach should provide a very powerful and accurate tool to predict the electrostatics of not only solution but also enzymatic reactions, as well as the solvation free energies of even larger systems, such as nucleic acid bases incorporated into DNA.

Analysis of catalytic mechanism of serine proteases. Viability of the ring-flip hypothesis

Quantum calculations are applied to the active site of serine proteases, including four specific residues and a water molecule, as well as a substrate and proton donors in the oxyanion hole. Residues are tethered to the protein backbone of an X-ray structure but otherwise allowed to move freely to their lowest energy positions. The viability of the ring-flip hypothesis, which proposes that a 180 degrees rotation of the His-57 imidazole ring facilitates the catalysis, is assessed by comparison of energies of configurations both before and after such a flip. Specifically considered is the contribution to catalysis of the Ser-214 residue and a water molecule that is observed in the active site. The calculations provide detailed information concerning the nature, geometry, and strength of hydrogen bonds that are formed within the active site at each stage of the enzymatic mechanism.

Substrate polarization in enzyme catalysis: QM/MM analysis of the effect of oxaloacetate polarization on acetyl-CoA enolization in citrate synthase

Citrate synthase is an archetypal carbon-carbon bond forming enzyme. It promotes the conversion of oxaloacetate (OAA) to citrate by catalyzing the deprotonation (enolization) of acetyl-CoA, followed by nucleophilic attack of the enolate form of this substrate on OAA to form a citryl-CoA intermediate and subsequent hydrolysis. OAA is strongly bound to the active site and its alpha-carbonyl group is polarized. This polarization has been demonstrated spectroscopically, [(Kurz et al., Biochemistry 1985;24:452-457; Kurz and Drysdale, Biochemistry 1987;26:2623-2627)] and has been suggested to be an important catalytic strategy. Substrate polarization is believed to be important in many enzymes. The first step, formation of the acetyl-CoA enolate intermediate, is thought to be rate-limiting in the mesophilic (pig/chicken) enzyme. We have examined the effects of substrate polarization on this key step using quantum mechanical/molecular mechanical (QM/MM) methods. Free energy profiles have been calculated by AM1/CHARMM27 umbrella sampling molecular dynamics (MD) simulations, together with potential energy profiles. To study the influence of OAA polarization, profiles were calculated with different polarization of the OAA alpha-carbonyl group. The results indicate that OAA polarization influences catalysis only marginally but has a larger effect on intermediate stabilization. Different levels of treatment of OAA are compared (MM or QM), and its polarization in the protein and in water analyzed at the B3LYP/6-31+G(d)/CHARMM27 level. Analysis of stabilization by individual residues shows that the enzyme mainly stabilizes the enolate intermediate (not the transition state) through electrostatic (including hydrogen bond) interactions: these contribute much more than polarization of OAA.

Ab initio QM/MM modelling of acetyl-CoA deprotonation in the enzyme citrate synthase

The first step of the reaction catalysed by the enzyme citrate synthase is studied here with high level combined quantum mechanical/molecular mechanical (QM/MM) methods (up to the MP2/6-31+G(d)//6-31G(d)/CHARMM level). In the first step of the reaction, acetyl-CoA is deprotonated by Asp375, producing an intermediate, which is the nucleophile for attack on the second substrate, oxaloacetate, prior to hydrolysis of the thioester bond of acetyl-CoA and release of the products. A central question has been whether the nucleophilic intermediate is the enolate of acetyl-CoA, the enol, or an 'enolic' intermediate stabilized by a 'low-barrier' hydrogen bond with His274 at the active site. The imidazole sidechain of His274 is neutral, and donates a hydrogen bond to the carbonyl oxygen of acetyl-CoA in substrate complexes. We have investigated the identity of the nucleophilic intermediate by QM/MM calculations on the substrate (keto), enolate, enol and enolic forms of acetyl-CoA at the active site of citrate synthase. The transition states for proton abstraction from acetyl-CoA by Asp375, and for transfer of the hydrogen bonded proton between His274 and acetyl-CoA have been modelled approximately. The effects of electron correlation are included by MP2/6-31G(d) and MP2/6-31+G(d) calculations on active site geometries produced by QM/MM energy minimization. The results do not support the hypothesis that a low-barrier hydrogen bond is involved in catalysis in citrate synthase, in agreement with earlier calculations. The acetyl-CoA enolate is identified as the only intermediate consistent with the experimental barrier for condensation, stabilized by conventional hydrogen bonds from His274 and a water molecule.

Theoretical studies of the transition-state structures and free energy barriers for base-catalyzed hydrolysis of amides

The transition-state structures and free energy barriers for the rate-determining step (i.e. the formation of a tetrahedral intermediate) of base-catalyzed hydrolysis of a series of amides in aqueous solution have been studied by performing first-principle electronic structure calculations using a hybrid supermolecule-polarizable continuum approach. The calculated results and a revisit of recently reported experimental proton inventory data reveal that the favorable transition-state structure optimized for the tetrahedral intermediate formation of hydroxide ion-catalyzed hydrolysis of formamide may have three solvating water molecules remaining on the attacking hydroxide oxygen and two additional water molecules attached to the carbonyl oxygen of formamide. The calculated results have also demonstrated interesting substituent effects on the optimized transition-state geometries, on the transition-state stabilization, and on the calculated free energy barriers for the base-catalyzed hydrolysis of amides. When some or all of the hydrogen atoms of formamide are replaced by methyl groups, the total number of water molecules hydrogen-bonding with the attacking hydroxide in the transition state decreases from three for formamide to two for N-methylacetamide, N,N-dimethylformamide (DMF), and N,N-dimethylacetamide (DMA). The larger substituents of the amide hinder the solvent water molecules approaching the attacking hydroxide oxygen in the transition state and, therefore, destabilize the transition-state structure and increase the free energy barrier. By using the optimized most favorable transition-state structures, the calculated free energy barriers, i.e., 21.6 (or 21.7), 22.7, 23.1, and 26.0 kcal/mol for formamide, N-methylacetamide, DMF, and DMA, respectively, are in good agreement with the available experimental free energy barriers, i.e., 21.2, 21.5, 22.6, and 24.1 kcal/mol for formamide, N-methylacetamide, DMF, and DMA, respectively.

Reaction path potential for complex systems derived from combined ab initio quantum mechanical and molecular mechanical calculations

Combined ab initio quantum mechanical and molecular mechanical calculations have been widely used for modeling chemical reactions in complex systems such as enzymes, with most applications being based on the determination of a minimum energy path connecting the reactant through the transition state to the product in the enzyme environment. However, statistical mechanics sampling and reaction dynamics calculations with a combined ab initio quantum mechanical (QM) and molecular mechanical (MM) potential are still not feasible because of the computational costs associated mainly with the ab initio quantum mechanical calculations for the QM subsystem. To address this issue, a reaction path potential energy surface is developed here for statistical mechanics and dynamics simulation of chemical reactions in enzymes and other complex systems. The reaction path potential follows the ideas from the reaction path Hamiltonian of Miller, Handy and Adams for gas phase chemical reactions but is designed specifically for large systems that are described with combined ab initio quantum mechanical and molecular mechanical methods. The reaction path potential is an analytical energy expression of the combined quantum mechanical and molecular mechanical potential energy along the minimum energy path. An expansion around the minimum energy path is made in both the nuclear and the electronic degrees of freedom for the QM subsystem internal energy, while the energy of the subsystem described with MM remains unchanged from that in the combined quantum mechanical and molecular mechanical expression and the electrostatic interaction between the QM and MM subsystems is described as the interaction of the MM charges with the QM charges. The QM charges are polarizable in response to the changes in both the MM and the QM degrees of freedom through a new response kernel developed in the present work. The input data for constructing the reaction path potential are energies, vibrational frequencies, and electron density response properties of the QM subsystem along the minimum energy path, all of which can be obtained from the combined quantum mechanical and molecular mechanical calculations. Once constructed, it costs much less for its evaluation. Thus, the reaction path potential provides a potential energy surface for rigorous statistical mechanics and reaction dynamics calculations of complex systems. As an example, the method is applied to the statistical mechanical calculations for the potential of mean force of the chemical reaction in triosephosphate isomerase.

The barrier for proton transport in aquaporins as a challenge for electrostatic models: The role of protein relaxation in mutational calculations

The origin of the barrier for proton transport through the aquaporin channel is a problem of general interest. It is becoming increasingly clear that this barrier is not attributable to the orientation of the water molecules across the channel but rather to the electrostatic penalty for moving the proton charge to the center of the channel. However, the reason for the high electrostatic barrier is still rather controversial. It has been argued by some workers that the barrier is due to the so-called NPA motif and/or to the helix macrodipole or to other specific elements. However, our works indicated that the main reason for the high barrier is the loss of the generalized solvation upon moving the proton charge from the bulk to the center of the channel and that this does not reflect a specific repulsive electrostatic interaction but the absence of sufficient electrostatic stabilization. At this stage it seems that the elucidation and clarification of the origin of the electrostatic barrier can serve as an instructive test case for electrostatic models. Thus, we reexamine the free-energy surface for proton transport in aquaporins using the microscopic free-energy perturbation/umbrella sampling (FEP/US) and the empirical valence bond/umbrella sampling (EVB/US) methods as well as the semimacroscopic protein dipole Langevin dipole model in its linear response approximation version (the PDLD/S-LRA). These extensive studies help to clarify the nature of the barrier and to establish the "reduced solvation effect" as the primary source of this barrier. That is, it is found that the barrier is associated with the loss of the generalized solvation energy (which includes of course all electrostatic effects) upon moving the proton charge from the bulk solvent to the center of the channel. It is also demonstrated that the residues in the NPA region and the helix dipole cannot be considered as the main reasons for the electrostatic barrier. Furthermore, our microscopic and semimacroscopic studies clarify the problems with incomplete alternative calculations, illustrating that the effects of various electrostatic elements are drastically overestimated by macroscopic calculations that use a low dielectric constant and do not consider the protein reorganization. Similarly, it is pointed out that microscopic potential of mean force calculations that do not evaluate the electrostatic barrier relative to the bulk water cannot be used to establish the origin of the electrostatic barrier. The relationship between the present study and calculations of pK(a)s in protein interiors is clarified, pointing out that approaches that are applied to study the aquaporin barrier should be validated by pK(a)s calculations. Such calculations also help to clarify the crucial role of solvation energies in establishing the barrier in aquaporins.

A self-stabilized model of the chymotrypsin catalytic pocket. The energy profile of the overall catalytic cycle

A model of the catalytic triad of chymotrypsin is built assuring the arrangement and properties as they are within the complete enzyme. The model contains 18 amino acid residues of chymotrypsin and its substrate. A total of 135 atoms (including 70 heavy atoms) were subjected to full ab initio geometry optimizations through 127 individual steps along the reaction coordinate of the complete catalytic mechanism. It was shown that the described model of the catalytic apparatus forms a self-stabilized molecule ensemble without the rest of the enzyme and substrate. According to the calculations, the formations of the first and second tetrahedral intermediates in the model have 20.3 and 15.7 kcal/mol activation energy barriers, respectively. Removing elements of the catalytic apparatus such as the (1) catalytic aspartate or (2) the anion hole, as well as (3) inserting a water molecule "early" in the catalytic process, or (4) introducing conformational rigidity of the substrate, results in an increase of the above energy barrier of the first catalytic step in the model by 6.4, 13.7, 3.7, and 4.1 kcal/mol, respectively. Based on the calculated process one can conclude that the catalytic reaction in this model is much more similar to the reaction in the enzyme than to the reference reaction. To our knowledge, this is the first model system that mimics the complete catalytic mechanism.

Probing the Structural and Electronic Factors Affecting the Adsorption and Reactivity of Alkenes in Acidic Zeolites Using DFT Calculations and Multivariate Statistical Methods

Quantum mechanical (QM) cluster calculations have been performed on a model of ZSM-5 at DFT and MP2 levels. We investigated how the adsorption energies and the energetics of alkoxide intermediate formation of six different alkene substrates, ethene, propene, 1-butene, cis/trans butene, and isobutene, vary in this zeolite model. An analysis of the DFT geometric, electronic, and energetic parameters of the zeolite-substrate complexes, transition states, and alkoxide intermediates is performed using principal components analysis (PCA) and partial least squares (PLS). These deliver an insight into the correlated changes that occur between molecular structure and energy along the reaction coordinate between the physisorbed and chemisorbed species within the zeolite. To the best of our knowledge, this is the first occasion multivariate techniques such as PCA or PLS have been employed to profile the changes in electronics, distances, and angles in QM calculations of catalytic systems such as zeolites. We find the calculated adsorption and the alkoxide intermediate energies correlate strongly with the absolute charge on the substrate and the length of the substrate double bond. The transition states' energies are not affected by the zeolite framework as modeled, which explains why they correlate strongly with the gas-phase substrate protonation energy. Our cluster results show that for ethene, propene, 1-butene, and isobutene, the relative energetics associated with the formation of the alkoxide intermediate in ZSM-5 follow the same trends as calculations where the effects of the framework are included.

Theoretical Studies on the Mechanism of Acid-Promoted Hydrolysis of N-Formylaziridine in Comparison with Formamide

[reaction: see text] We present an ab initio study of the acid-promoted hydrolysis reaction mechanism of N-formylaziridine in comparison with formamide. Since the rate of amide hydrolysis reactions depends on the formation of the tetrahedral intermediate, we focused our attention mainly on the reactant complex, the tetrahedral intermediate, and the transition state connecting these two stationary points. Geometries were optimized using the density functional theory, and the energetics were refined using ab initio theory including electron correlation. Solvent effects were investigated by using polarizable continuum method calculations. The proton-transfer reaction between the O-protonated and N-protonated amides was investigated. In acidic media, despite that the N-protonated species is more stable than the O-protonated one, it is predicted that both N-protonated and O-protonated pathways compete in the hydrolysis reaction of N-formylaziridine.

Theoretical analysis of peptidyl alpha-ketoheterocyclic inhibitors of human neutrophil elastase: Insight into the mechanism of inhibition and the application of QM/MM calculations in structure-based drug design

It has been suggested from QSAR data (P. D. Edwards, D. J. Wolanin, D.A. Andisik and M. W. Davis, J. Med. Chem., 1995, 38, 76) that the inhibition of elastase by peptidyl alpha-ketoheterocyclic inhibitors can occur in two ways, the less potent inhibitors forming a non-bonded Michaelis complex and the more potent set a covalently bonded enzyme-substrate intermediate. We report QM/MM studies of both binding and reactivity that confirm these findings, showing that the activity of the least potent set of inhibitors correlates with the calculated binding energy, and that of the more potent set correlates with the stability of the intermediate. These calculations show that QM/MM methods can be successfully employed to understand complicated structure-activity relationships and might be employed in the design and assessment of new inhibitors.

Computer simulations of enzyme catalysis: methods, progress, and insights

Understanding the action of enzymes on an atomistic level is one of the important aims of modern biophysics. This review describes the state of the art in addressing this challenge by simulating enzymatic reactions. It considers different modeling methods including the empirical valence bond (EVB) and more standard molecular orbital quantum mechanics/molecular mechanics (QM/MM) methods. The importance of proper configurational averaging of QM/MM energies is emphasized, pointing out that at present such averages are performed most effectively by the EVB method. It is clarified that all properly conducted simulation studies have identified electrostatic preorganization effects as the source of enzyme catalysis. It is argued that the ability to simulate enzymatic reactions also provides the chance to examine the importance of nonelectrostatic contributions and the validity of the corresponding proposals. In fact, simulation studies have indicated that prominent proposals such as desolvation, steric strain, near attack conformation, entropy traps, and coherent dynamics do not account for a major part of the catalytic power of enzymes. Finally, it is pointed out that although some of the issues are likely to remain controversial for some time, computer modeling approaches can provide a powerful tool for understanding enzyme catalysis.

Modulation of the NO trans effect in heme proteins: implications for the activation of soluble guanylate cyclase

The binding of NO to the iron heme in guanylate cyclase and other heme proteins induces the cleavage of the proximal histidine bonded to the metal. In this study we assess by means of density functional theory (DFT) electronic structure calculations the role of H-bonding to histidine in the modulation of this effect. We have considered in the first place a model of the isolated active site coordinated with imidazole and imidazolate to mimic the effects of a very strong H-bond. We have also investigated four selected ferrous heme proteins with different proximal histidine environments: the O(2) sensing FixL, horseradish peroxidase C, and the alpha and beta subunits of human hemoglobin. Our results indicate that polarization and charge transfer effects associated with H-bonding to the proximal histidine play a fundamental role in the modulation of the NO trans effect in heme proteins. We also find computational evidence suggesting that protein structural constraints may affect significantly the cleavage of the Fe-His bond.

The intramolecular mechanism for the second proton transfer in triosephosphate isomerase (TIM): a QM/FE approach

The intramolecular mechanism we earlier proposed [Alagona, G.; Desmeules, P.; Ghio, C.; Kollman, P. A. J Am Chem Soc 1984, 106, 3623] for the second proton transfer of the reaction catalyzed by triosephosphate isomerase (TIM) is examined ab initio at the HF and MP2/6-31+G** levels in vacuo for two conformers of the enediolate phosphate (ENEP), with the ethereal oxygen of the phosphate group either syn (X), as in the crystal structure, or anti (Y) with respect to the enediolate carbonyl O. The barrier height for the intramolecular proton transfer occurring in enediolate is very sensitive to electron correlation corrections. The MP2 internal energy barrier is much lower than the HF one, while the free energy (FE) barrier is even more favorable, indicating that the enzyme presence is not requested to speed up that step. An investigation of the dynamical aspects of the mechanism, along the pathway from ENEP A (with H on O(1)) to TS and from TS to ENEP B (with H on O(2)), was, however, carried out in the presence of the enzyme field while using a neutral His-95 with its proton on Ndelta. To perform the FE simulations, it was necessary to parametrize in the AMBER force-field the ENEP A, TS and B species, whose partial charges have been determined with the RESP procedure, with the X and Y arrangements of the phosphate head. Actually, the FE/QM approach produced a low barrier and a substantial balance between A and B, especially at the MP2 level. The trajectories, analyzed paying a particular attention to the positions assumed by His-95 and by the other active site residues, put forward somewhat different H-bond patterns around the X or Y enediolate phosphate.

Evolution of catalytic antibody repertoire in autoimmune mice

We have attempted to efficiently obtain catalytic antibodies (catAbs) with amidase/esterase activity in the expanded sequence space of the antibody repertoire. In doing so, we used an autoimmune mouse strain, MRL/lpr, that is known to produce enhanced levels of autoantibodies. We applied different types of haptens, such as, and, that are supposed to mimic the transition state of the substrate in the ester/amide hydrolysis. Among them, hapten (2) could not be used, as it was readily broken down after synthesis. Upon immunization with hapten (1), catAbs preferentially evolved in MRL/lpr mice, but this did not happen upon immunization with haptens (3) and (4). Independently, immunization to MRL/lpr mice with successfully elicited the catAbs with the ability to activate vitamin B(6) prodrugs. The common observation seen in these two cases is that most of the catAbs derived from MRL/lpr mice by hapten (1) and half of them by hapten (5) had a Lys at H95, which is at the junctional N region between the V(H) and J(H) gene segments. Despite the conservation of Lys (H95), analyses of the N-region and utilization of the D gene segment in the heavy chain gene showed that these catAbs were from several independent clones of the same family. Studies of site-directed mutagenesis suggest that, in the catAbs elicited from hapten (1), a Lys (H95) and a His (L91) are involved in the catalytic function. Both residues are known to interact with the phosphonate moiety of hapten (1). Such studies also suggest that, in the catAbs elicited from hapten (5), a Lys (H95) and a His (H35) are involved in the catalytic function. These basic amino acids seem to be important for binding to the phosphonate hapten, as they were not changed even after extensive evolution following multiple mutations. By contrast, in normal BALB/c mice, immunization of hapten (1) resulted in eliciting catAbs in lower yield and the majority were the non-catAbs, whose sequences were quite different from those of the catAbs from MRL/lpr mice. They were clonally related to one another and most of them originated from a single clone. The positions of the interacting key residues in the CDRs that interact with the phosphorus moiety strongly differ between our catAbs and other reported catAbs with esterase/amidase activity, which were elicited by the phosphonate/phosphonamidate haptens from normal mice. Further comparison of antibodies elicited by the phosphorus haptens, such as DNA, RNA, phosphocholine, and phosphotyrosine, indicated that none of them had sequence similarity in the basic amino acids and their positions in the CDRs, except for one example, which is anti-DNA antibody elicited from C3H-lpr mice. Analysis based on the classification of canonical structures of the antibodies again suggested that our catAbs derived from MRL/lpr mice belong to an unusual class that is not listed in the literature. Taken together, the above evidence suggests that the unique catalytic subsets that existed in the initial repertoire in the MRL/lpr mice could effectively be captured by the phosphonate haptens through the interaction with the Lys at H95. In the BALB/c mice, however, another noncatalytic subset with an ability to bind only to a moiety other than the phosphonate moiety alternatively evolved, because of the lowest abundance or elimination of the catalytic subsets.

Quantum mechanical methods for enzyme kinetics

This review discusses methods for the incorporation of quantum mechanical effects into enzyme kinetics simulations in which the enzyme is an explicit part of the model. We emphasize three aspects: (a) use of quantum mechanical electronic structure methods such as molecular orbital theory and density functional theory, usually in conjunction with molecular mechanics; (b) treating vibrational motions quantum mechanically, either in an instantaneous harmonic approximation, or by path integrals, or by a three-dimensional wave function coupled to classical nuclear motion; (c) incorporation of multidimensional tunneling approximations into reaction rate calculations.

Modeling enzyme reactivity in organic solvents and water through computer simulations

In this article, we review how molecular modeling techniques can be used to shed light on how water and organic solvents influence the reactivity of enzymes. The application of thermodynamics-based models allowed the first qualitative predictions on the selectivity of many reaction types. However, it was with the application of quantum mechanical/molecular mechanical (QM/MM) methods that quantitative models of actual reactivity patterns could be realistically formulated.