Density functional based parametrization of a valence bond method and its applications in quantum-classical molecular dynamics simulations of enzymatic reactions

An approximate valence bond (AVB) method was parametrized at a microscopic level for proton transfer and hydroxyanion nucleophilic reactions in enzyme catalytic processes. The method was applied to describe hydrolytic activity of phospholipase A₂. The AVB parametrization is based on density functional and conventional ab initio calculations calibrated with respect to experimental data in the gas phase. The method was used as a fast generator of the potential energy function in a quantum-classical molecular dynamics (QCMD) simulations describing atomic motions as well as propagation of the proton wave function in the enzyme active site. The protein environment surrounding the active site and solvent effects are included in the model via electrostatic interactions perturbing the original AVB Hamiltonian. © 1996 John Wiley & Sons, Inc.     

Nuclear quantum effects on an enzyme-catalyzed reaction with reaction path potential: proton transfer in triosephosphate isomerase

Nuclear quantum mechanical effects have been examined for the proton transfer reaction catalyzed by triosephosphate isomerase, with the normal mode centroid path integral molecular dynamics based on the potential energy surface from the recently developed reaction path potential method. In the simulation, the primary and secondary hydrogens and the C and O atoms involving bond forming and bond breaking were treated quantum mechanically, while all other atoms were dealt classical mechanically. The quantum mechanical activation free energy and the primary kinetic isotope effects were examined. Because of the quantum mechanical effects in the proton transfer, the activation free energy was reduced by 2.3 kcal/mol in comparison with the classical one, which accelerates the rate of proton transfer by a factor of 47.5. The primary kinetic isotope effects of kH/kD and kH/kT were estimated to be 4.65 and 9.97, respectively, which are in agreement with the experimental value of 4+/-0.3 and 9. The corresponding Swain-Schadd exponent was predicted to be 3.01, less than the semiclassical limit value of 3.34, indicating that the quantum mechanical effects mainly arise from quantum vibrational motion rather than tunneling. The reaction path potential, in conjunction with the normal mode centroid molecular dynamics, is shown to be an efficient computational tool for investigating the quantum effects on enzymatic reactions involving proton transfer.     

Fragmentation of conjugated amides at the C?C(O) bond in electrospray mass spectrometry: a proton‐bound dimeric intermediate identified by the kinetic method

Fragmentation of protonated amides in mass spectrometry at the C? C(O) bond, which competes with the C(O)? N bond cleavage, was observed when the amide group is conjugated with an unsaturated moiety. In the case of N‐methylbenzamides bearing different substituents, this reaction gives rise to protonated methyl isocyanate and protonated benzenes. The kinetic method was applied to this reaction in spite of the fact that this is a fragmentation of a molecular species (rather than proton‐bound dimers). A correlation is established between the intensities of the two product ions and the proton affinities (PAs) of the corresponding fragment molecules, which is similar to that when the kinetic method is used in the determination of PAs on proton‐bound dimers. This result provides strong evidence that the reaction proceeds via a proton‐bound methyl isocyanate/benzene complex intermediate. In addition, the PA of isocyanic acid, which is involved in the fragmentation of unsubstituted benzamide, is compared with that of benzene and a downward revised value, 175 kcal/mol, is recommended. Copyright © 2004 John Wiley & Sons, Ltd.     

Ab initio study of β-lactam antibiotics. I. Potential energy surface for the amidic CN bond breaking in the β-lactam + OH− reaction

The potential energy surface of the β-lactam + OH^? reaction, related to the mode of action of β-lactam antibiotics, was investigated using the ab initio Hartree—Fock method with the STO-3G basis set. Three possible reaction paths for the B_AC2 breaking of the amidic CN bond were obtained and discussed. The minimum-energy reaction path is characterized by the following processes: (1) the formation of a tetrahedral intermediate, ≈ 121 kcal mol^?1 more stable than the reagents; (2) a barrier, ≈ 15 kcal mol^?1 above the intermediate, which is mainly due to the partial breaking of the amidic bond; (3) the complete breaking of the amidic bond concerted with a proton transfer till the formation of the final product, ≈ 34 kcal mol^?1 more stable than the intermediate. The evolution of some molecular orbitals and of the electron population along the reaction path was also discussed.     

Excited States of Conjugated Hydrocarbons Using the Molecular Mechanics - Valence Bond (MMVB) Method: Conical Intersections and Dynamics

We developed the molecular mechanics—valence bond (MMVB) method following an original suggestion of Jean-Paul Malrieu and coworkers. By coupling a parameterized Heisenberg Hamiltonian to a standard classical force field (MM2), reliable ground and excited state geometries of conjugated hydrocarbons can be rapidly optimized. The MMVB method was central to our development of algorithms for locating conical intersections and calculating their associated decay dynamics. Here, we briefly review the chemical applications of MMVB to date, and present two new studies using the photostability of pyracylene and the excited state decay dynamics of the photochromic dihydroazulene/vinylheptafulvene (DHA/VHF) reaction.     

Potential energy functions for an intramolecular proton transfer reaction in the ground and excited state

We describe the development of empirical potential functions for the study of the excited state intramolecular proton transfer reaction in 1-(trifuloroacetylamino)-naphtaquinone (TFNQ). The potential is a combination of the standard CHARMM27 force field for the backbone structure of TFNQ and an empirical valence bond formalism for the proton transfer reaction. The latter is parameterized to reproduce the potential energies both in the ground and the excited state, determined at the CASPT2 level of theory. Parameters describing intermolecular interactions are fitted to reproduce molecular dipole moments computed at the CASSCF level of theory and to reproduce ab initio hydrogen bonding energies and geometries for TFNQ-water bimolecular complexes. The utility of this potential energy function was examined by computing the potentials of mean force for the proton transfer reactions in the gas phase and in water, in both electronic states. The ground state PMF exhibits little solvent effects, whereas computed potential of mean force shows a solvent stabilization of 2.5 kcal mol⁻¹ in the product state region, suggesting proton transfer is more pronounced in polar solvents, consistent with experimental findings. Electronic Supplementary Material The online version of this article (doi:) contains supplementary material, which is available to authorized users. Contribution to the Fernando Bernardi Memorial Issue.     

Computational analysis of the proton translocation from Asp96 to schiff base in bacteriorhodopsin

The potential energy change during the M --> N process in bacteriorhodopsin has been evaluated by ab initio quantum chemical and advanced quantum chemical calculations following molecular dynamics (MD) simulations. Many previous experimental studies have suggested that the proton transfer from Asp96 to the Schiff base occurs under the following two conditions: (1) the hydrogen bond between Thr46 and Asp96 breaks and Thr46 is detached from Asp96 and (2) a stable chain of four water molecules spans an area from Asp96 --> Schiff base. In this work, we successfully reproduced the proton-transfer process occurring under these two conditions by molecular dynamics and quantum chemical calculations. The quantum chemical computation revealed that the proton transfer from Asp96 to Shiff base occurs in two-step reactions via an intermediate in which an H(3)O(+) appears around Ala215. The activation energy for the proton transfer in the first reaction was calculated to be 9.7 kcal/mol, which enables fast and efficient proton pump action. Further QM/MM (quantum mechanical/molecular mechanical) and FMO (fragment molecular orbital) calculations revealed that the potential energy change during the proton transfer is tightly regulated by the composition and the geometry of the surrounding amino acid residues of bacteriorhodopsin. Here, we report in detail the Asp96 --> Schiff base proton translocation mechanism of bacteriorhodopsin. Additionally, we discuss the effectiveness of combining quantum chemical calculations with truncated cluster models followed by advanced quantum chemical calculations applied to a whole protein to elucidate its reaction mechanism.     

Theoretical study on the mechanism of a ring-opening reaction of oxirane by the active-site aspartic dyad of HIV-1 protease

Two possible mechanisms of the irreversible inhibition of HIV-1 protease by epoxide inhibitors are investigated on an enzymatic model using ab initio (MP2) and density functional theory (DFT) methods (B3LYP, MPW1K and M05-2X). The calculations predict the inhibition as a general acid-catalyzed nucleophilic substitution reaction proceeding by a concerted SN2 mechanism with a reaction barrier of ca. 15-21 kcal mol(-1). The irreversible nature of the inhibition is characterized by a large negative reaction energy of ca. -17-(-24) kcal mol(-1). A mechanism with a direct proton transfer from an aspartic acid residue of the active site onto the epoxide ring has been shown to be preferred compared to one with the proton transfer from the acid catalyst facilitated by a bridging catalytic water molecule. Based on the geometry of the transition state, structural data important for the design of irreversible epoxide inhibitors of HIV-1 protease were defined. Here we also briefly discuss differences between the epoxide ring-opening reaction in HIV-1 protease and epoxide hydrolase, and the accuracy of the DFT method used.     

A theoretical study of carbon-carbon bond formation by a Michael-type addition

A theoretical study of the Michael-type addition of 1,3-dicarbonyl compounds to α,β-unsaturated carbonyl compounds has been performed in the gas phase by means of the AM1 semiempirical method and by density functional theory (DFT) calculations within the B3LYP and M06-2X hybrid functionals. A molecular model has been selected to mimic the role of a base, which is traditionally used as a catalyst in Michael reactions, an acetate moiety to modulate its basicity, and point charges to imitate the stabilization of the negative charge developed in the substrate during the reaction when taking place in enzymatic environments. Results of the study of six different reactions obtained at the three different levels of calculations show that the reaction takes place in three steps: in the first step the α proton of the acetylacetone is abstracted by the base, then the nucleophilic attack on the β-carbon of the α,β-unsaturated carbonyl compound takes place generating the negatively charged enolate intermediate, and finally the product is formed through a proton transfer back from the protonated base. According to the energy profiles, the rate limiting step corresponds to the abstraction of the proton or the carbon-carbon bond formation step, depending on substituents of the substrates and method of calculation. The effect of the substituents on the acidity of the α proton of the acetylacetone and the steric hindrance can be analyzed by comparing these two separated steps. Moreover, the result of adding a positive charge close to the center that develops a negative charge during the reaction confirms the catalytic role of the oxyanion hole proposed in enzyme catalysed Michael-type additions. Stabilization of the intermediate implies, in agreement with the Hammond postulate, a reduction of the barrier of the carbon-carbon bond formation step. Our results can be used to predict the features that a new designed biocatalyst must present to efficiently accelerate this fundamental reaction in organic synthesis.     

A QM/MM study of the racemization of vinylglycolate catalyzed by mandelate racemase enzyme.

The experimentally postulated mechanism for the interconversion between (S)-vinylglycolate and (R)-vinylglycolate catalyzed by mandelate racemase enzyme consists of a two-step quite symmetric process through a dianionic enolic intermediate that is formed after the abstraction of the alpha-proton of vinylglycolate by a basic enzymatic residue and is then reprotonated by another residue. The challenging problem behind this reaction is how the enzyme manages to stabilize such an intermediate, that is, how it lowers enough the high pK(a) of the alpha-proton for the reaction to take place. The QM/MM simulations performed in this paper indicate that catalysis is based on the stabilization of the negative charge developed on the substrate along the reaction. We have identified three different reaction mechanisms starting from different quasi-degenerate structures of the substrate-enzyme complex. In two of them the stabilizing role is done by means of a catalytic proton transfer that avoids the formation of a dianionic intermediate, and they involve six steps instead of the two experimentally proposed. On the contrary, the third mechanism passes through a dianionic species stabilized by the concerted approach of a protonated enzymatic residue during the proton abstraction. The potential energy barriers theoretically found along these mechanisms are qualitatively in good agreement with the experimental free energy barriers determined for racemization of vinylglycolate and mandelate. The theoretical study of the effect of the mutation of Glu317 by Gln317 in the kinetics of the reaction reveals the important role in the catalysis of the hydrogen bond formed by Glu317 in the native enzyme, as only one of the mechanisms, the slower one, is able to produce the racemization in the active site of the mutant. However, we have found that this hydrogen bond is not an LBHB within our model.     

Ab initio study of β-lactam antibiotics. II. Potential energy surface for the amidic CN bond breaking in the 3-cephem + OH− reaction and comparison with the β-lactam + OH− reaction

The potential energy surface for the β-lactam amidic CN bond breaking in the 3-cephem + OH^? reaction was investigated by using the ab initio Hartree—Fock method with a 9s6p/7s3p/3s basis set. The investigated reaction is a model of the reaction between an antibiotic cephalosporin and an enzymatic nucleophilic group, this last reaction being related to the mode of action of β-lactam antibiotics. The minimum-energy reaction path is characterized by a tetrahedral intermediate ≈ 116 kcal mol^?1 more stable than the reagents, by a barrier which corresponds to the partial breaking of the amidic bond and is ≈ 7 kcal mol^?1 above the intermediate and by a product ≈ 31 kcal mol^?1 more stable than the intermediate. The analysis of the wavefunction along the reaction path and the comparison with the β-lactam + OH^? reaction pointed out the role of electron-withdrawing groups on the height of the barrier and the role of intramolecular hydrogen bonds on the structure and energy of the product. The calculations suggest a model of the antibiotic activity of cephalosporins which is compared with previous qualitative pictures.     

Investigation of the mechanism of the cell wall DD-carboxypeptidase reaction of penicillin-binding protein 5 of Escherichia coli by quantum mechanics/molecular mechanics calculations

Penicillin-binding protein 5 (PBP 5) of Escherichia coli hydrolyzes the terminal D-Ala-D-Ala peptide bond of the stem peptides of the cell wall peptidoglycan. The mechanism of PBP 5 catalysis of amide bond hydrolysis is initial acylation of an active site serine by the peptide substrate, followed by hydrolytic deacylation of this acyl-enzyme intermediate to complete the turnover. The microscopic events of both the acylation and deacylation half-reactions have not been studied. This absence is addressed here by the use of explicit-solvent molecular dynamics simulations and ONIOM quantum mechanics/molecular mechanics (QM/MM) calculations. The potential-energy surface for the acylation reaction, based on MP2/6-31+G(d) calculations, reveals that Lys47 acts as the general base for proton abstraction from Ser44 in the serine acylation step. A discrete potential-energy minimum for the tetrahedral species is not found. The absence of such a minimum implies a conformational change in the transition state, concomitant with serine addition to the amide carbonyl, so as to enable the nitrogen atom of the scissile bond to accept the proton that is necessary for progression to the acyl-enzyme intermediate. Molecular dynamics simulations indicate that transiently protonated Lys47 is the proton donor in tetrahedral intermediate collapse to the acyl-enzyme species. Two pathways for this proton transfer are observed. One is the direct migration of a proton from Lys47. The second pathway is proton transfer via an intermediary water molecule. Although the energy barriers for the two pathways are similar, more conformers sample the latter pathway. The same water molecule that mediates the Lys47 proton transfer to the nitrogen of the departing D-Ala is well positioned, with respect to the Lys47 amine, to act as the hydrolytic water in the deacylation step. Deacylation occurs with the formation of a tetrahedral intermediate over a 24 kcal x mol(-1) barrier. This barrier is approximately 2 kcal x mol(-1) greater than the barrier (22 kcal x mol(-1)) for the formation of the tetrahedral species in acylation. The potential-energy surface for the collapse of the deacylation tetrahedral species gives a 24 kcal x mol(-1) higher energy species for the product, signifying that the complex would readily reorganize and pave the way for the expulsion of the product of the reaction from the active site and the regeneration of the catalyst. These computational data dovetail with the knowledge on the reaction from experimental approaches.     

Catalytic role of metal ion in the selection of competing reaction paths: a first principles molecular dynamics study of the enzymatic reaction in ribozyme

By using finite temperature first principles molecular dynamics, the mechanism of the enzymatic reaction of ribozyme was investigated for both the anionic and the radical charge states of the modeled RNA fragment. In the case of the anionic system, a pseudorotation and the subsequent 3' --> 2' migration occur in a vacuum, rather than the self-cleavage of the phosphodiester. On the other hand, when either a divalent metal ion (Mg(2+)) catalyst or the continuous hydrogen bond network of the solvent is present, the reaction path of the anionic species changes dramatically, going toward the transesterification channel. In a radical system, the transesterification can occur without a metal catalyst, as a consequence of the displacement of a hole (empty electronic state) along the reaction path. Thus, the present analysis suggests that a metal ion might be essential not only in lowering the activation barrier but also in selecting the reaction path among those corresponding to possible different charge states of the intermediate structure in vivo. Furthermore, simulation of the anionic species in solution shows that, in the absence of a metal catalyst, water molecules cooperate with the proton transfer via a proton wire mechanism and the hydrogen bond network plays a crucial role in preventing pseudorotations. On the other hand, when a metal cation is present in the vicinity of the site where the nucleophilic attack occurs, the hydrogen bond network is interrupted and detachment of the proton, enhanced by the catalyst, does not give rise to any proton-transfer process.     

Enzymatic assay technique for the determination of aspartame

An enzymatic assay technique was developed for the determination of the artificial sweetener aspartame. The peptide bond of aspartame was first cleaved by peptidase to release aspartic acid. In the presence of α-ketoglutarate, aspartic acid was then transaminated by aspartate aminotransferase to glutamate. The reaction was monitored by following the oxygen consumption during the enzymatic oxidation of glutamate by glutamate oxidase. A linear relationship between oxygen consumption and aspartame concentration up to 200 μM was obtained. The assay technique was applicable to the determination of aspartame in a variety of dietary food products. The results obtained agreed well with those determined by liquid chromatography and those reported by the product manufacturers.     

A computational study of the role of hydrogen bonds in SN1 and E1 reactions

The reaction between tertiary butyl chloride and water clusters was examined by applying density functional theory calculations. The carbonium ion t-Bu(+) that is normally sandwiched between the water clusters was found to be absent, such that a Cbond;O covalent bond was formed in the intermediate (Int1) after heterolysis. An (H(2)O)(4) cluster is able to bridge the front and rear of the central carbon and promotes heterolysis. A correlation between bond interchanges at the central carbon and proton relays is presented. Stereochemical scrambling in the solvolysis products is discussed in terms of this correlation. In addition, an E1 pathway for the elimination product, iso-butene, is found from Int1.     

Car-Parrinello and path integral molecular dynamics study of the hydrogen bond in the chloroacetic acid dimer system

We have studied the double proton transfer (DPT) reaction in the cyclic dimer of chloroacetic acid using both classical and path integral Car-Parrinello molecular dynamics. We also attempt to quantify the errors in the potential energy surface that arise from the use of a pure density functional. In the classical dynamics a clear reaction mechanism can be identified, where asynchronized DPT arises due to coupling between the O-H stretching oscillator and several low energy intermolecular vibrational modes. This mechanism is considerably altered when quantum tunneling is permitted in the simulation. The introduction of path integrals leads to considerable changes in the thermally averaged molecular geometry, leading to shorter and more centered hydrogen bond linkages.     

Approach to nonadiabatic transitions by density matrix evolution and molecular dynamics simulations

Many biological processes are characterized by an essentially quantum dynamical event, such as electron or proton transfer, in a complex classical environment. To treat such processes properly by computer simulation, allowing nonadiabatic transitions involving excited states, we recently developed a density matrix evolution (DME ) method [H. J. C. Berendsen and J. Mavri, J. Phys. Chem, 97 , 13464 (1993)] which simulates the dynamics of quantum systems embedded in a classical environment. The formalism of the method is presented and an overview of the applications ranging from collisions of a quantum harmonic oscillator with noble gas atoms to proton tunneling in a double-well hydrogen bond is given. The methodology for treatment of proton-transfer processes with inclusion of excited states is presented. Future application of the method on biologically interesting processes, such as proton transfer in enzymatic reactions, is discussed. © 1996 John Wiley & Sons, Inc.     

Promoting vibrations in human purine nucleoside phosphorylase. A molecular dynamics and hybrid quantum mechanical/molecular mechanical study

Crystallographic studies of human purine nucleoside phosphorylase (hPNP) with several transition-state (TS) analogues in the immucillin family showed an unusual geometric arrangement of the atoms O-5', O-4', and O(P), the nucleophilic phosphate oxygen, lying in a close three-oxygen stack. These observations were corroborated by extensive experimental kinetic isotope effect analysis. We propose that protein-facilitated dynamic modes in hPNP cause this stack, centered on the ribosyl O-4' oxygen, to squeeze together and push electrons toward the purine ring, stabilizing the oxacarbenium character of the TS. As the N-ribosidic bond is cleaved during the reaction, the pK(a) values of N-7 and O-6 increase by the electron density expelled by the oxygen-stack compression toward the purine ring. Increased electron density in the purine ring improves electrostatic interactions with nearby residues and facilitates the abstraction of a proton from a solvent proton or an unidentified general acid, making the purine a better leaving group, and accelerating catalysis. Classical and mixed quantum/classical molecular dynamics (MD) simulations of the Michaelis complex of hPNP with the substrates guanosine and phosphate were performed to assess the existence of protein-promoting vibrations (PPVs). Analogous simulations were performed for the substrates in aqueous solution. In the catalytic site, the O-5', O-4', and O(P) oxygens vibrate at frequencies of ca. 125 and 465 cm(-1), as opposed to 285 cm(-1) in the absence of hPNP. The hybrid quantum mechanical/molecular mechanical method was used to assess whether this enzymatic vibration pushing the oxygens together is coupled to the reaction coordinate, and thus has a direct positive impact on catalysis. The potential energy surface for the phosphorolysis reaction for several snapshots taken from the classical MD simulation showed substantial differences in oxygen compression. Our calculations showed the existence of PPVs coupled to the reaction coordinate, which effect electronic alterations in the active site by pushing the three oxygen centers together in proximity, and accelerate substrate turnover in the phosphorolysis reaction catalyzed by hPNP.     

Catalytic mechanism of guanine deaminase: an ONIOM and molecular dynamics study

The catalytic mechanism of Bacillus subtilis guanine deaminase (bGD), a Zn metalloenzyme, has been investigated by a combination of quantum mechanical calculations using the multilayered ONIOM method and molecular dynamics simulations. In contrast to a previously proposed catalytic mechanism, which requires the bound guanine to assume a rare tautomeric state, the ONIOM calculations showed that the active-site residues of the enzyme do not affect the tautomeric state of guanine, and consequently the bound guanine is a tautomer that is the most abundant in aqueous solution. Two residues, Glutamate 55 and Aspartate 114, were found to play important roles in proton shuttling in the reaction. The proposed reaction path is initiated by proton transfer from a Zn-bound water to protonate Asp114. This process may be quite complex and rather dynamic in nature, as revealed by the molecular dynamics (MD) simulations, whereby another water may bridge the Zn-bound water and Asp114, which then is eliminated by positioning of guanine in the active site. The binding of guanine stabilizes protonated Asp114 by hydrogen bond formation. Asp114 can then transfer its proton to the N3 of the bound guanine, facilitating the nucleophilic attack on C2 of the guanine by the Zn-bound hydroxide to form a tetrahedral intermediate. This occurs with a rather low barrier. Glu55 then transfers a proton from the Zn-hydroxide to the amino group of the reaction intermediate and, at this point, the C2-N2 bond has lengthened by 0.2 A compared to guanine, making C2-N2 bond cleavage more facile. The C2-N2 bond breaks forming ammonia, with an energy barrier of approximately 8.8 kcal/mol. Ammonia leaves the active site, and xanthine is freed by the cleavage of the Zn-O2 bond, with a barrier approximately 8.4 kcal/mol. Along this reaction path, the highest barrier comes from C2-N2 bond cleavage, while the barrier from the cleavage of the Zn-O2 bond is slightly smaller. The Zn-O2 bond can be broken without the assistance of water during the release of xanthine.     

SINDO1 study of the photoreaction of tetramethylene sulfone

The mechanism of the photoreaction of tetramethylene sulfone (TMSO₂) was investigated by the semiempirical molecular orbital (MO) method SINDO1. The relevant low-lying potential energy surfaces, which were calculated with limited configuration interaction (CI), were studied by optimizing intermediates and transition structures and by introducing linear interpolations between these stationary points. The main initial reaction step for all important products is an α cleavage of one C? S bond. This leads to an intermediate that can be classified as an excited singlet diradical. Its electronic structure is described with a two-electron, three-orbital model. Starting from this initially generated intermediate, the reaction branches into several pathways leading to various products. Feasible reaction pathways were established for all important products consistent with experiments. © 1995 by John Wiley & Sons, Inc.