A comparative study of trypsin specificity based on QM/MM molecular dynamics simulation and QM/MM GBSA calculation

Hydrogen bonding and polar interactions play a key role in identification of protein-inhibitor binding specificity. Quantum mechanics/molecular mechanics molecular dynamics (QM/MM MD) simulations combined with DFT and semi-empirical Hamiltonian (AM1d, RM1, PM3, and PM6) methods were performed to study the hydrogen bonding and polar interactions of two inhibitors BEN and BEN1 with trypsin. The results show that the accuracy of treating the hydrogen bonding and polar interactions using QM/MM MD simulation of PM6 can reach the one obtained by the DFT QM/MM MD simulation. Quantum mechanics/molecular mechanics generalized Born surface area (QM/MM-GBSA) method was applied to calculate binding affinities of inhibitors to trypsin and the results suggest that the accuracy of binding affinity prediction can be significantly affected by the accurate treatment of the hydrogen bonding and polar interactions. In addition, the calculated results also reveal the binding specificity of trypsin: (1) the amidinium groups of two inhibitors generate favorable salt bridge interaction with Asp189 and form hydrogen bonding interactions with Ser190 and Gly214, (2) the phenyl of inhibitors can produce favorable van der Waals interactions with the residues His58, Cys191, Gln192, Trp211, Gly212, and Cys215. This systematic and comparative study can provide guidance for the choice of QM/MM MD methods and the designs of new potent inhibitors targeting trypsin.     

QM/MM free energy simulations: recent progress and challenges

Due to the higher computational cost relative to pure molecular mechanical (MM) simulations, hybrid quantum mechanical/molecular mechanical (QM/MM) free energy simulations particularly require a careful consideration of balancing computational cost and accuracy. Here, we review several recent developments in free energy methods most relevant to QM/MM simulations and discuss several topics motivated by these developments using simple but informative examples that involve processes in water. For chemical reactions, we highlight the value of invoking enhanced sampling technique (e.g. replica-exchange) in umbrella sampling calculations and the value of including collective environmental variables (e.g. hydration level) in metadynamics simulations; we also illustrate the sensitivity of string calculations, especially free energy along the path, to various parameters in the computation. Alchemical free energy simulations with a specific thermodynamic cycle are used to probe the effect of including the first solvation shell into the QM region when computing solvation free energies. For cases where high-level QM/MM potential functions are needed, we analyse two different approaches: the QM/MM–MFEP method of Yang and co-workers and perturbative correction to low-level QM/MM free energy results. For the examples analysed here, both approaches seem productive although care needs to be exercised when analysing the perturbative corrections.     

Pentacoordinated phosphorus revisited by high‐level QM/MM calculations

Enzymes catalyzing phosphoryl transfer reactions are extremely efficient and are involved in crucial biochemical processes. The mechanisms of these enzymes are complex due to the diversity of substrates that are involved. The reaction can proceed through a pentacoordinated phosphorus species that is either a stable intermediate or a transition state (TS). Because of this, the first X‐ray structure of a pentacoordinated phosphorus intermediate in the β‐phosphoglucomutase enzyme aroused great interest but also much controversy. To provide new insights into the nature of that structure, we have determined the reaction path of the phosphorylation step using high‐level QM/MM calculations, and have also calculated the geometry of a complex with a transition state analogue (TSA) that has been suggested to be the actual species in the crystal. The protein crystalline environment has been modeled so as to mimic the experimental conditions. We conclude that the pentacoordinated phosphorus formed in this enzyme is not a stable species but a TS, which gives an activation energy for phosphorylation in agreement with kinetic results. We also show that the TSA is a good mimic of the true TS. We have performed a new crystallographic refinement of the original diffraction map of the pentacoordinated phosphorus structure with the MgF TSA. The new fit improves significantly with respect to the original one, which strongly supports that Allen and coworkers wrongly assigned the X‐ray structure to a pentavalent phosphorane. Proteins 2010. © 2010 Wiley‐Liss, Inc.     

A comparative study of galactose oxidase and active site analogs based on QM/MM Car-Parrinello simulations

A parallel study of the radical copper enzyme galactose oxidase (GOase) and a low molecular weight analog of the active site was performed with dynamical density functional and mixed quantum-classical calculations. This combined approach enables a direct comparison of the properties of the biomimetic and the natural systems throughout the course of the catalytic reaction. In both cases, five essential forms of the catalytic cycle have been investigated: the resting state in its semi-reduced (catalytically inactive) and its oxidized (catalytically active) form, A ^semi and A ^ox, respectively; a protonated intermediate B; the transition state for the rate-determining hydrogen abstraction step C, and its product D. For A and B the electronic properties of the biomimetic compound are qualitatively very similar to the ones of the natural target. However, in agreement with the experimentally observed difference in catalytic activity, the calculated activation energy for the hydrogen abstraction step is distinctly lower for GOase (16 kcal/mol) than for the mimetic compound (21 kcal/mol). The enzymatic transition state is stabilized by a delocalization of the unpaired spin density over the sulfur-modified equatorial tyrosine Tyr272, an effect that for geometric reasons is essentially absent in the biomimetic compound. Further differences between the mimic and its natural target concern the structure of the product of the abstraction step, which is characterized by a weakly coordinated aldehyde complex for the latter and a tightly bound linear complex for the former. Received 14 October 1999 · Accepted: 19 January 2000     

A QM/MM study of the complexes formed by aluminum and iron with serum transferrin at neutral and acidic pH

Serum transferrin (sTf) transports iron in serum and internalizes in cells via receptor mediated endocytosis. Additionally, sTf has been identified as the predominant aluminum carrier in serum. Some questions remain unclear about the exact mechanism for the metal release or whether the aluminum and iron show the same binding mode during the entire process. In the present work, simulation techniques at quantum and atomic levels have been employed in order to gain access into a molecular level understanding of the metal-bound sTf complex, and to describe the binding of Al(III) and Fe(III) ions to sTf. First, hybrid quantum mechanics/molecular mechanics (QM/MM) molecular dynamics (MD) simulations were carried out in order to analyze the dynamics of the aluminum-loaded complex, taking into account the different pH conditions in blood and into the cell. Moreover, the complexes formed by transferrin with Al(III) and Fe(III) were optimized with high level density functional theory (DFT)/MM methods. All these results indicate that the interaction mode of Al(III) and Fe(III) with sTf change upon different pH conditions, and that the coordination of Al(III) and Fe(III) is not equivalent during the metal intake, transport and release processes. Our results emphasize the importance of the pH on the metal binding and release mechanism and suggest that Al(III) can follow the iron pathway to get access into cells, although once there, it may show a different binding mode, leading to a different mechanism for its release.     

QM/MM investigation of structure and spectroscopic properties of a vanadium-containing peroxidase

We present a comprehensive analysis of the most likely ground state configuration of the resting state of vanadium dependent chloroperoxidase (VCPO) based on quantum mechanics/molecular mechanics (QM/MM) evaluations of ground state properties, UV-vis spectra and NMR chemical shifts. Within the QM/MM framework, density functional theory (DFT) calculations are used to characterize the resting state of VCPO via time-dependent density functional theory (TD-DFT) calculations of electronic excitation energies and NMR chemical shifts. Comparison with available experimental data allows us to determine the most likely protonation state of VCPO, a state which results in a doubly protonated axial oxygen, a site largely stabilized by hydrogen bonds. We found that the bulk of the protein that is beyond the immediate layer surrounding the cofactor, has an important electrostatic effect on the absorption maximum. Through examination of frontier orbitals, we analyze the nature of two bound water molecules and the extent to which relevant residues in the active site influence the spectroscopy calculations.     

Reaction mechanism of caspases: insights from QM/MM Car-Parrinello simulations

Caspases are fundamental targets for pharmaceutical interventions in a variety of diseases involving disregulated apoptosis. Here, we present a quantum mechanics/molecular mechanics Car-Parrinello study of key steps of the enzymatic reaction for a representative member of this family, caspase-3. The hydrolysis of the acyl-enzyme complex is described at the density functional (BLYP) level of theory while the protein frame and solvent are treated using the GROMOS96 force field. These calculations show that the attack of the hydrolytic water molecule implies an activation free energy of ca. DeltaF(A) approximately equal 19 +/- 4 kcal/mol in good agreement with experimental data and leads to a previously unrecognized gem-diol intermediate that can readily (DeltaF(A) approximately equal 5 +/- 3 kcal/mol) evolve to the enzyme products. Our findings assist in elucidating the striking difference in catalytic activity between caspases and other structurally well-characterized cysteine proteases (papains and cathepsins) and may help design novel transition-state analog inhibitors.     

A QM/MM study of the catalytic mechanism of SAM methyltransferase RlmN from Escherichia coli

RlmN is a radical S‐adenosylmethionine (SAM) enzyme that catalyzes the C2 methylation of adenosine 2503 (A2503) in 23S rRNA and adenosine 37 (A37) in several Escherichia coli transfer RNAs (tRNA). The catalytic reaction of RlmN is distinctly different from that of typical SAM‐dependent methyltransferases that employs an S_N2 mechanism, but follows a ping‐pong mechanism which involves the intermediate methylation of a conserved cysteine residue. Recently, the x‐ray structure of a key intermediate in the RlmN reaction has been reported, allowing us to perform combined quantum mechanics and molecular mechanics (QM/MM) calculations to delineate the reaction details of RlmN at atomic level. Starting from the Cross‐Linked RlmN C118A?tRNA complex, the possible mechanisms for both the formation and the resolution of the cross‐linked species (IM2) have been illuminated. On the basis of our calculations, IM2 is formed by the attack of the C355‐based methylene radical on the sp²‐hybridized C2 of the adenosine ring, corresponding to energy barrier of 14.4 kcal/mol, and the resolution of IM2 is confirmed to follow a radical fragmentation mechanism. The cleavage of C′–S′ bond of mC355‐A37 cross‐link is in concert with the deprotonation of C2 by C118 residue, which is the rate‐limiting step with an energy barrier of 17.4 kcal/mol. Moreover, the cleavage of C′–S′ bond of IM2 can occur independently, that is, it does not require the loss of an electron of IM2 and the formation of disulfide bond between C355 and C118 as precondition. These findings would deepen the understanding of the catalysis of RlmN.     

The QM/MM molecular dynamics and free energy simulations of the acylation reaction catalyzed by the serine-carboxyl peptidase kumamolisin-As

Quantum mechanical/molecular mechanical molecular dynamics and free energy simulations are performed to study the acylation reaction catalyzed by kumamolisin-As, a serine-carboxyl peptidase, and to elucidate the catalytic mechanism and the origin of substrate specificity. It is demonstrated that the nucleophilic attack by the serine residue on the substrate may not be the rate-limiting step for the acylation of the GPH*FF substrate. The present study also confirms the earlier suggestions that Asp164 acts as a general acid during the catalysis and that the electrostatic oxyanion hole interactions may not be sufficient to lead a stable tetrahedral intermediate along the reaction pathway. Moreover, Asp164 is found to act as a general base during the formation of the acyl-enzyme from the tetrahedral intermediate. The role of dynamic substrate assisted catalysis (DSAC) involving His at the P1 site of the substrate is examined for the acylation reaction. It is demonstrated that the bond-breaking and -making events at each stage of the reaction trigger a change of the position for the His side chain and lead to the formation of the alternative hydrogen bonds. The back and forth movements of the His side chain between the C=O group of Pro at P2 and Odelta2 of Asp164 in a ping-pong-like mechanism and the formation of the alternative hydrogen bonds effectively lower the free energy barriers for both the nucleophilic attack and the acyl-enzyme formation and may therefore contribute to the relatively high activity of kumamolisin-As toward the substrates with His at the P1 site.     

QM/MM study of D-fructose in aqueous solution

The QM/MM molecular dynamics methodology was applied to the study of the two main D-fructose tautomers present in aqueous solution, beta-D-fructofuranose and beta-D-fructopyranose. The solute was treated at the AM1 semi-empirical level, and for the solvent water molecules we used the TIP3P potential. We analyzed the structure of the water molecules around the hydroxyl groups to explain the differences in sweet taste between the two tautomers.     

Asymmetric biodegradation of the nerve agents Sarin and VX by human dUTPase: chemometrics,molecular docking and hybrid QM/MM calculations

Organophosphorus compounds (OP) nerve agents are among the most toxic chemical substances known. Their toxicity is due to their ability to bind to acetylcholinesterase. Currently, some enzymes, such as phosphotriesterase, human serum paraoxonase 1 and diisopropyl fluorophosphatase, capable of degrading OP, have been characterized. Regarding the importance of bioremediation methods for detoxication of OP, this work aims to study the interaction modes between the human human deoxyuridine triphosphate nucleotidohydrolase (dUTPase) and Sarin and VX, considering their R_p and S_p enantiomers, to evaluate the asymmetric catalysis of those compounds. In previous work, this enzyme has shown good potential to degrade phosphotriesters, and based on this characteristic, we have applied the human dUTPase to the OP degradation. Molecular docking, chemometrics and mixed quantum and molecular mechanics calculations have been employed, showing a good interaction between dUTPase and OP. Two possible reaction mechanisms were tested, and according to our theoretical results, the catalytic degradation of OP by dUTPase can take place via both mechanisms, beyond being stereoselective, that is, dUTPase cleaves one enantiomer preferentially in relation to other. Chemometric techniques provided excellent assistance for performing this theoretical investigation. The dUTPase study shows importance by the fact of it being a human enzyme.

Communicated by Ramaswamy H. Sarma     

Mechanisms of guanosine triphosphate hydrolysis by Ras and Ras-GAP proteins as rationalized by ab initio QM/MM simulations

The hydrolysis reaction of guanosine triphosphate (GTP) by p21(ras) (Ras) has been modeled by using the ab initio type quantum mechanical-molecular mechanical simulations. Initial geometry configurations have been prompted by atomic coordinates of the crystal structure (PDBID: 1QRA) corresponding to the prehydrolysis state of Ras in complex with GTP. Multiple searches of minimum energy geometry configurations consistent with the hydrogen bond networks have been performed, resulting in a series of stationary points on the potential energy surface for reaction intermediates and transition states. It is shown that the minimum energy reaction path is consistent with an assumption of a two-step mechanism of GTP hydrolysis. At the first stage, a unified action of the nearest residues of Ras and the nearest water molecules results in a substantial spatial separation of the gamma-phosphate group of GTP from the rest of the molecule (GDP). This phase of hydrolysis process proceeds through the low barrier (16.7 kcal/mol) transition state TS1. At the second stage, the inorganic phosphate is formed in consequence of proton transfers mediated by two water molecules and assisted by the Gln61 residue from Ras. The highest transition state at this segment, TS3, is estimated to have an energy 7.5 kcal/mol above the enzyme-substrate complex. The results of simulations are compared to the previous findings for the GTP hydrolysis in the Ras-GAP (p21(ras)-p120(GAP)) protein complex. Conclusions of the modeling lead to a better understanding of the anticatalytic effect of cancer causing mutation of Gln61 from Ras, which has been debated in recent years.     

Role of the arginine finger in Ras.RasGAP revealed by QM/MM calculations

In the Ras.Ras.GAP complex, hydrolysis of guanosine triphosphate is strongly accelerated GAP as compared to Ras alone. This is largely attributed to the arginine finger R789(GAP) pointing to AlF(x) in the transition state analogue. We performed QM/MM simulations where triphosphate was treated using the quantum mechanical method of density functional theory, while the protein complex and water environment were described classically using MD. Compared to Ras, the crucial electron shift, bond stretching and distortion towards an eclipsed gamma-to-beta orientation are much more pronounced. The arginine finger is shown to act by displacing water out of the binding niche. The resulting enhanced electrostatic field catalyses the cleavage step.     

Thermodynamics calculation of protein–ligand interactions by QM/MM polarizable charge parameters

The calculation of protein–ligand binding free energy (ΔG) is of great importance for virtual screening and drug design. Molecular dynamics (MD) simulation has been an attractive tool to investigate this scientific problem. However, the reliability of such approach is affected by many factors including electrostatic interaction calculation. Here, we present a practical protocol using quantum mechanics/molecular mechanics (QM/MM) calculations to generate polarizable QM protein charge (QMPC). The calculated QMPC of some atoms in binding pockets was obviously different from that calculated by AMBER ff03, which might significantly affect the calculated ΔG. To evaluate the effect, the MD simulations and MM/GBSA calculation with QMPC for 10 protein–ligand complexes, and the simulation results were then compared to those with the AMBER ff03 force field and experimental results. The correlation coefficient between the calculated ΔΔG using MM/GBSA under QMPC and the experimental data is .92, while that with AMBER ff03 force field is .47 for the complexes formed by streptavidin or its mutants and biotin. Moreover, the calculated ΔΔG with QMPC for the complexes formed by ERβ and five ligands is positively related to experimental result with correlation coefficient of .61, while that with AMBER ff03 charge is negatively related to experimental data with correlation coefficient of .42. The detailed analysis shows that the electrostatic polarization introduced by QMPC affects the electrostatic contribution to the binding affinity and thus, leads to better correlation with experimental data. Therefore, this approach should be useful to virtual screening and drug design.     

QM/MM studies on the catalytic mechanism of Phenylethanolamine N-methyltransferase

Epinephrine is a naturally occurring adrenomedullary hormone that transduces environmental stressors into cardiovascular actions. As the only route in the catecholamine biosynthetic pathway, Phenylethanolamine N-methyltransferase (PNMT) catalyzes the synthesis of epinephrine. To elucidate the detailed mechanism of enzymatic catalysis of PNMT, combined quantum-mechanical/molecular-mechanical (QM/MM) calculations were performed. The calculation results reveal that this catalysis contains three elementary steps: the deprotonation of protonated norepinphrine, the methyl transferring step and deprotonation of the methylated norepinphrine. The methyl transferring step was proved to be the rate-determining step undergoing a SN2 mechanism with an energy barrier of 16.4 kcal/mol. During the whole catalysis, two glutamic acids Glu185 and Glu219 were proved to be loaded with different effects according to the calculations results of the mutants. These calculation results can be used to explain the experimental observations and make a good complementarity for the previous QM study.     

QM/MM study of the mechanism of enzymatic limonene 1,2-epoxide hydrolysis

Limonene 1,2-epoxide hydrolase (LEH) is completely different from those of classic epoxide hydrolases (EHs) which catalyze the hydrolysis of epoxides to vicinal diols. A novel concerted general acid catalysis step involving the Asp101-Arg99-Asp132 triad is proposed to play an important role in the mechanism. Combined quantum-mechanical/molecular-mechanical (QM/MM) calculations gave activation barriers of 16.9 and 25.1 kcal/mol at the B3LYP/6-31G(d,p)//CHARMM level for nucleophilic attack on the more and less substituted epoxide carbons, respectively. Furthermore, the important roles of residues Arg99, Tyr53 and Asn55 on mutated LEH were evaluated by QM/MM-scanned energy mapping. These results may provide an explanation for site-directed mutagenesis.     

QM/MM study of the C—C coupling reaction mechanism of CYP121, an essential Cytochrome p450 of Mycobacterium tuberculosis

Among 20 p450s of Mycobacterium tuberculosis (Mt), CYP121 has received an outstanding interest, not only due to its essentiality for bacterial viability but also because it catalyzes an unusual carbon–carbon coupling reaction. Based on the structure of the substrate bound enzyme, several reaction mechanisms were proposed involving first Tyr radical formation, second Tyr radical formation, and C—C coupling. Key and unknown features, being the nature of the species that generate the first and second radicals, and the role played by the protein scaffold each step. In the present work we have used classical and quantum based computer simulation methods to study in detail its reaction mechanism. Our results show that substrate binding promotes formation of the initial oxy complex, Compound I is the responsible for first Tyr radical formation, and that the second Tyr radical is formed subsequently, through a PCET reaction, promoted by the presence of key residue Arg386. The final C—C coupling reaction possibly occurs in bulk solution, thus yielding the product in one oxygen reduction cycle. Our results thus contribute to a better comprehension of MtCYP121 reaction mechanism, with direct implications for inhibitor design, and also contribute to our general understanding of these type of enzymes. Proteins 2014; 82:1004–1021. © 2013 Wiley Periodicals, Inc.     

QM/MM modeling the Ras-GAP catalyzed hydrolysis of guanosine triphosphate

The mechanism of the hydrolysis reaction of guanosine triphosphate (GTP) by the protein complex Ras-GAP (p21(ras) - p120(GAP)) has been modeled by the quantum mechanical-molecular mechanical (QM/MM) and ab initio quantum calculations. Initial geometry configurations have been prompted by atomic coordinates of a structural analog (PDBID:1WQ1). It is shown that the minimum energy reaction path is consistent with an assumption of two-step chemical transformations. At the first stage, a unified motion of Arg789 of GAP, Gln61, Thr35 of Ras, and the lytic water molecule results in a substantial spatial separation of the gamma-phosphate group of GTP from the rest of the molecule (GDP). This phase of hydrolysis process proceeds through the low-barrier transition state TS1. At the second stage, Gln61 abstracts and releases protons within the subsystem including Gln61, the lytic water molecule and the gamma-phosphate group of GTP through the corresponding transition state TS2. Direct quantum calculations show that, in this particular environment, the reaction GTP + H(2)O --> GDP + H(2)PO(4) (-) can proceed with reasonable activation barriers of less than 15 kcal/mol at every stage. This conclusion leads to a better understanding of the anticatalytic effect of cancer-causing mutations of Ras, which has been debated in recent years.