Molecular mechanisms of antibiotic resistance: QM/MM modelling of deacylation in a class A beta-lactamase

Modelling of the first step of the deacylation reaction of benzylpenicillin in the E. coli TEM1 beta-lactamase (with B3LYP/6-31G + (d)//AM1-CHARMM22 quantum mechanics/molecular mechanics methods) shows that a mechanism in which Glu166 acts as the base to deprotonate a conserved water molecule is both energetically and structurally consistent with experimental data; the results may assist the design of new antibiotics and beta-lactamase inhibitors.     

Ab initio QM/MM study of class A beta-lactamase acylation: dual participation of Glu166 and Lys73 in a concerted base promotion of Ser70

Beta-lactamase acquisition is the most prevalent basis for Gram-negative bacteria resistance to the beta-lactam antibiotics. The mechanism used by the most common class A Gram-negative beta-lactamases is serine acylation followed by hydrolytic deacylation, destroying the beta-lactam. The ab initio quantum mechanical/molecular mechanical (QM/MM) calculations, augmented by extensive molecular dynamics simulations reported herein, describe the serine acylation mechanism for the class A TEM-1 beta-lactamase with penicillanic acid as substrate. Potential energy surfaces (based on approximately 350 MP2/6-31+G calculations) reveal the proton movements that govern Ser70 tetrahedral formation and then collapse to the acyl-enzyme. A remarkable duality of mechanism for tetrahedral formation is implicated. Following substrate binding, the pathway initiates by a low energy barrier (5 kcal mol(-1)) and an energetically favorable transfer of a proton from Lys73 to Glu166, through the catalytic water molecule and Ser70. This gives unprotonated Lys73 and protonated Glu166. Tetrahedral formation ensues in a concerted general base process, with Lys73 promoting Ser70 addition to the beta-lactam carbonyl. Moreover, the three-dimensional potential energy surface also shows that the previously proposed pathway, involving Glu166 as the general base promoting Ser70 through a conserved water molecule, exists in competition with the Lys73 process. The existence of two routes to the tetrahedral species is fully consistent with experimental data for mutant variants of the TEM beta-lactamase.     

Theoretical modeling of the heme group with a hybrid QM/MM method

The quality of the results obtained in calculations with the hybrid QM/MM method IMOMM on systems where the heme group is partitioned in QM and MM regions is evaluated through the performance of calculations on the 4‐coordinate [Fe(P)] (P = porphyrin), the 5‐coordinate [Fe(P)(1−(Me)Im)] (Im = imidazole) and the 6‐coordinate [Fe(P)(1−(Me)Im)(O₂)] systems. The results are compared with those obtained from much more expensive pure quantum mechanics calculations on model systems. Three different properties are analyzed—namely, the optimized geometries, the binding energies of the axial ligands to the heme group, and the energy cost of the biochemically relevant out‐of‐plane displacement of the iron atom. Agreement is especially good in the case of optimized geometries and energy cost of out‐of‐plane displacements, with larger discrepancies in the case of binding energies. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 282–294, 2000     

The reaction mechanism of hydroxyethylphosphonate dioxygenase: a QM/MM study

By employing ab initio quantum mechanical/molecular mechanical (QM/MM) and molecular dynamics (MD) simulations, we have provided further evidence against the previously proposed hydroperoxylation or hydroxylation mechanism of hydroxyethylphosphonate dioxygenase (HEPD). HEPD employs an interesting catalytic cycle based on concatenated bifurcations. The first bifurcation is based on the abstraction of hydrogen atoms from the substrate, which leads to a distal or proximal hydroperoxo species (Fe-OOH or Fe-(OH)O). The second and the third bifurcations refer to the carbon-carbon bond cleavage reaction. And this is achieved through a tridentate intermediate, or employing a proton-shuttle assisted mechanism, in which the residue Glu(176) or the Fe(IV)=O group serves as a general base. The reaction directions seem to be tunable and show significant environment dependence. This mechanism can provide a comprehensive interpretation for the seemingly contradicting experimental evidences and provide insight into the development of biochemistry and material sciences.     

Antibiotic deactivation by a dizinc beta-lactamase: mechanistic insights from QM/MM and DFT studies

Hybrid quantum mechanical/molecular mechanical (QM/MM) methods and density functional theory (DFT) were used to investigate the initial ring-opening step in the hydrolysis of moxalactam catalyzed by the dizinc L1 beta-lactamase from Stenotrophomonas maltophilia. Anchored at the enzyme active site via direct metal binding as suggested by a recent X-ray structure of an enzyme-product complex (Spencer, J.; et al. J. Am. Chem. Soc. 2005, 127, 14439), the substrate is well aligned with the nucleophilic hydroxide that bridges the two zinc ions. Both QM/MM and DFT results indicate that the addition of the hydroxide nucleophile to the carbonyl carbon in the substrate lactam ring leads to a metastable intermediate via a dominant nucleophilic addition barrier. The potential of mean force obtained by SCC-DFTB/MM simulations and corrected by DFT/MM calculations yields a reaction free energy barrier of 23.5 kcal/mol, in reasonable agreement with the experimental value of 18.5 kcal/mol derived from kcat of 0.15 s(-1). It is further shown that zinc-bound Asp120 plays an important role in aligning the nucleophile, but accepts the hydroxide proton only after the nucleophilic addition. The two zinc ions are found to participate intimately in the catalysis, consistent with the proposed mechanism. In particular, the Zn(1) ion is likely to serve as an "oxyanion hole" in stabilizing the carbonyl oxygen, while the Zn(2) ion acts as an electrophilic catalyst to stabilize the anionic nitrogen leaving group.     

Mixed quantum mechanical/molecular mechanical (QM/MM) study of the deacylation reaction in a penicillin binding protein (PBP) versus in a class C beta-lactamase

The origin of the substantial difference in deacylation rates for acyl-enzyme intermediates in penicillin-binding proteins (PBPs) and beta-lactamases has remained an unsolved puzzle whose solution is of great importance to understanding bacterial antibiotic resistance. In this work, accurate, large-scale mixed ab initio quantum mechanical/molecular mechanical (QM/MM) calculations have been used to study the hydrolysis of acyl-enzyme intermediates formed between cephalothin and the dd-peptidase of Streptomyces sp. R61, a PBP, and the Enterobacter cloacae P99 cephalosporinase, a class C beta-lactamase. Qualitative and, in the case of P99, quantitative agreement was achieved with experimental kinetics. The faster rate of deacylation in the beta-lactamase is attributed to a more favorable electrostatic environment around Tyr150 in P99 (as compared to that for Tyr159 in R61) which facilitates this residue's function as the general base. This is found to be in large part accomplished by the ability of P99 to covalently bind the ligand without concurrent elimination of hydrogen bonds to Tyr150, which proves not to be the case with Tyr159 in R61. This work provides an essential foundation for further work in this area, such as selecting mutations capable of converting the PBP into a beta-lactamase.     

Theoretical modeling of open-shell molecules in solution: a QM/MM molecular dynamics approach

In this work, the GLOB model, an effective and reliable computational approach well suited for ab initio and QM/MM molecular dynamics simulations of complex molecular systems in solution, has been applied to study two representative open-shell systems, the cobalt(II) ion and the glycine radical in aqueous solution, with special reference to their structural and magnetic properties. The main structural features of the solvent cage around the cobalt ion and the hydrogen bonding patterns around the neutral and zwitterionic forms of the glycine radical have been investigated in some detail. The general good agreement with experiments supports the use of the present model to investigate more challenging and biological/technological relevant open-shell systems.     

pKa calculations in solution and proteins with QM/MM free energy perturbation simulations: a quantitative test of QM/MM protocols

The accuracy of biological simulations depends, in large part, on the treatment of electrostatics. Due to the availability of accurate experimental values, calculation of pKa provides stringent evaluation of computational methods. The generalized solvent boundary potential (GSBP) and Ewald summation electrostatic treatments were recently implemented for combined quantum mechanical and molecular mechanics (QM/MM) simulations by our group. These approaches were tested by calculating pKa shifts due to differences in electronic structure and electrostatic environment; the shifts were determined for a series of small molecules in solution, using various electrostatic treatments, and two residues (His 31, Lys 102) in the M102K T4-lysozyme mutant with large pKa shifts, using the GSBP approach. The calculations utilized a free energy perturbation scheme with the QM/MM potential function involving the self-consistent charge density functional tight binding (SCC-DFTB) and CHARMM as the QM and MM methods, respectively. The study of small molecules demonstrated that inconsistent electrostatic models produced results that were difficult to correct in a robust manner; by contrast, extended electrostatics, GSBP, and Ewald simulations produced consistent results once a bulk solvation contribution was carefully chosen. In addition to the electrostatic treatment, the pKa shifts were also sensitive to the level of the QM method and the scheme of treating QM/MM Coulombic interactions; however, simple perturbative corrections based on SCC-DFTB/CHARMM trajectories and higher level single point energy calculations were found to give satisfactory results. Combining all factors gave a root-mean-square difference of 0.7 pKa units for the relative pKa values of the small molecules compared to experiment. For the residues in the lysozyme, an accurate pKa shift was obtained for His 31 with multiple nanosecond simulations. For Lys 102, however, the pKa shift was estimated to be too large, even after more than 10 nanosecond simulations for each lambda window; the difficulty was due to the significant, but slow, reorganization of the protein and water structure when Lys 102 was protonated. The simulations support that Lys 102 is deprotonated in the X-ray structure and the protein is highly destabilized when this residue is protonated.     

Binding affinities by alchemical perturbation using QM/MM with a large QM system and polarizable MM model

The most general way to improve the accuracy of binding‐affinity calculations for protein–ligand systems is to use quantum‐mechanical (QM) methods together with rigorous alchemical‐perturbation (AP) methods. We explore this approach by calculating the relative binding free energy of two synthetic disaccharides binding to galectin‐3 at a reasonably high QM level (dispersion‐corrected density functional theory with a triple‐zeta basis set) and with a sufficiently large QM system to include all short‐range interactions with the ligand (744–748 atoms). The rest of the protein is treated as a collection of atomic multipoles (up to quadrupoles) and polarizabilities. Several methods for evaluating the binding free energy from the 3600 QM calculations are investigated in terms of stability and accuracy. In particular, methods using QM calculations only at the endpoints of the transformation are compared with the recently proposed non‐Boltzmann Bennett acceptance ratio (NBB) method that uses QM calculations at several stages of the transformation. Unfortunately, none of the rigorous approaches give sufficient statistical precision. However, a novel approximate method, involving the direct use of QM energies in the Bennett acceptance ratio method, gives similar results as NBB but with better precision, ～3 kJ/mol. The statistical error can be further reduced by performing a greater number of QM calculations. © 2015 Wiley Periodicals, Inc.     

A QM/MM study of a nucleophilic aromatic substitution reaction catalyzed by 4-chlorobenzoyl-CoA dehalogenase

Calculated using a QM/MM method, the free energy profile for the conversion of 4-chlorobenzoate to 4-hydroxybenzoate catalyzed by 4-chlorobenzoyl-CoA dehalogenase indicates the existence of a stable Meisenheimer complex.     

Macrophomate synthase: QM/MM simulations address the Diels-Alder versus Michael-Aldol reaction mechanism

Macrophomate synthase (MPS) of the phytopathogenic fungus Macrophoma commelinae catalyzes the transformation of 2-pyrone derivatives to the corresponding benzoate analogues. The transformation proceeds through three separate chemical reactions, including decarboxylation of oxalacetate to produce pyruvate enolate, two C-C bond formations between 2-pyrone and pyruvate enolate that form a bicyclic intermediate, and final decarboxylation with concomitant dehydration. Although some evidence suggests that the second step of the reaction catalyzed by MPS is a Diels-Alder reaction, definite proof that the C-C bond formations occur via a concerted mechanism is still required. An alternative route for formation of the C-C bonds is a stepwise Michael-aldol reaction. In this work, mixed quantum and molecular mechanics (QM/MM) combined with Monte Carlo simulations and free-energy perturbation (FEP) calculations were performed to investigate the relative stabilities of the transition states (TS) for both reaction mechanisms. The key results are that the Diels-Alder TS model is 17.7 and 12.1 kcal/mol less stable than the Michael and aldol TSs in the stepwise route, respectively. Therefore, this work indicates that the Michael-aldol mechanism is the route used by MPS to catalyze the second step of the overall transformation, and that the enzyme is not a natural Diels-Alderase, as claimed by Ose and co-workers (Nature 2003, 422, 185-189; Acta Crystallogr. 2004, D60, 1187-1197). A modified link-atom treatment for the bonds at the QM/MM interface is also presented.     

A critical evaluation of different QM/MM frontier treatments with SCC-DFTB as the QM method

The performance of different link atom based frontier treatments in QM/MM simulations was evaluated critically with SCC-DFTB as the QM method. In addition to the analysis of gas-phase molecules as in previous studies, an important element of the present work is that chemical reactions in realistic enzyme systems were also examined. The schemes tested include all options available in the program CHARMM for SCC-DFTB/MM simulation, which treat electrostatic interactions due to the MM atoms close to the QM/MM boundary in different ways. In addition, a new approach, the divided frontier charge (DIV), has been implemented in which the partial charge associated with the frontier MM atom ("link host") is evenly distributed to the other MM atoms in the same group. The performance of these schemes was evaluated based on properties including proton affinities, deprotonation energies, dipole moments, and energetics of proton transfer reactions. Similar to previous work, it was found that calculated proton affinities and deprotonation energies of alcohols, carbonic acids, amino acids, and model DNA bases are very sensitive to the link atom scheme; the commonly used single link atom approach often gives error on the order of 15 to 20 kcal/mol. Other schemes give better and, on average, mutually comparable results. For proton transfer reactions, encouragingly, both activation barriers and reaction energies are fairly insensitive (within a typical range of 2-4 kcal/mol) to the link atom scheme due to error cancellation, and this was observed for both gas-phase and enzyme systems. Therefore, the effect of using different link atom schemes in QM/MM simulations is rather small for chemical reactions that conserve the total charge. Although the current study used an approximate DFT method as the QM level, the observed trends are expected to be applicable to QM/MM methods with use of other QM approaches. This observation does not mean to encourage QM/MM simulations without careful benchmark in the study of specific systems, rather it emphasizes that other technical details, such as the treatment of long-range electrostatics, tend to play a more important role and need to be handled carefully.     

Chaperoned alchemical free energy simulations: a general method for QM, MM, and QM/MM potentials

A general method for alchemical free energy simulations using QM, MM, and QM/MM potential is developed by introducing "chaperones" to restrain the structures, particularly near the end points. A calculation of the free energy difference between two triazole tautomers in aqueous solution is used to illustrate the method.     

Flexible-boundary QM/MM calculations: II. Partial charge transfer across the QM/MM boundary that passes through a covalent bond

Recently, based on the principle of electronic chemical potential equalization and the principle of charge conservation, we proposed a flexible-boundary scheme that allows both partial charge transfer and self-consistent polarization between the quantum mechanical (QM) and molecular mechanical (MM) subsystems in QM/MM calculations; the scheme was applied to study the atomic charges in selected ion–solvent complexes. In the present contribution, we further extend the flexible-boundary treatment to handle the QM/MM boundary passing through covalent bonds. We find that the flexible-boundary redistributed charge and dipole schemes yield reasonable agreement with full-QM calculations for a number of molecular ions and amino acids with charged side chains. Using the full-QM results as reference, the mean unsigned deviations are computed to be 0.06 e for atomic partial charges of the QM atoms, 0.11 e for the amounts of charge transfer between the QM and MM subsystems, and 0.016 Å for the lengths of the covalent bonds that directly connect the QM and MM subsystems. The results indicate the importance of accounting for partial charge transfer across the QM/MM boundary when the QM subsystems are charged.     

烟酰胺酶催化水解脱氨反应的QM/MM分子动力学模拟

采用QM(PM3)/MM分子动力学(MD)方法模拟了烟酰胺酶催化烟酰胺水解脱氨形成烟碱酸的反应过程.计算结果表明,硫的亲核进攻是整个反应的速率控制步骤.当改变Ala155所在Loop区的位置,在亲核进攻时,底物能够自由旋转,可以加速亲核进攻过程并降低整个催化反应的能垒.讨论了氨分子的离去过程和质子传递过程的相关细节.为烟酰胺酶的定点突变以及脱氨酶的设计提供了有益的参考.     

抗原肽与MHC分子相互作用QM/MM分子动力学模拟研究

在MHC I类(major histocompatibility complex class I)分子抗原加工提呈过程中抗原蛋白在抗原提呈细胞(antigen-presenting cells, APC)的胞浆中被蛋白酶体(proteasome)裂解成短肽peptide, 由转运相关蛋白(transporter associated with antigen processing, TAP)将蛋白酶体裂解产生的短肽片段从胞浆转运至内质网腔. 短肽peptide在内质网中与新生成的MHC I类分子结合, 形成peptide-MHC复合体被提呈到APC细胞表面, 与T细胞表面抗原受体(T cell receptor, TCR)特异性识别结合, 使得CTL细胞开始活化、增殖、分化, 进而对肿瘤细胞进行特异性杀伤. 目前对CTL细胞如何识别抗原肽-MHC复合物分子及抗原短肽peptide如何与主要组织相容性复合体MHC分子的相互作用识别结合的机理还不是很清楚. 传统的预测CTL细胞表位的方法没有考虑受体与配体结合过程中电子结构的变化, 电子结构的变化需要用量子力学方法来处理. 本文采用QM/MM多尺度生物大分子的分子动力学模拟方法, 以天然抗原肽TAX (LLFGYPVYVYU)与HLA-A*0201分子结合的晶体结构为模板, 替换抗原肽“锚点”氨基酸, 将口袋氨基酸残基的原子极化电荷在空间形成的静电势用电多极矩分量表示. 用箱线图分析每个口袋氨基酸分子静电势变化和功能, 确定Pocket B的Glu63和Lys66的功能是精细识别氨基酸和一级结合氨基酸, Pocket F的Asp77, Tyr84的功能是精细识别氨基酸, 而Asp77, Lys146是一级结合氨基酸, 表明QM/MM方法在提取抗原肽与MHC I类分子识别结合特异性信息是可行的, 这对了解免疫识别机理和指导肿瘤疫苗的开发都具有指导意义.     

QuanPol: A full spectrum and seamless QM/MM program

The quantum chemistry polarizable force field program (QuanPol) is implemented to perform combined quantum mechanical and molecular mechanical (QM/MM) calculations with induced dipole polarizable force fields and induced surface charge continuum solvation models. The QM methods include Hartree–Fock method, density functional theory method (DFT), generalized valence bond theory method, multiconfiguration self‐consistent field method, Møller–Plesset perturbation theory method, and time‐dependent DFT method. The induced dipoles of the MM atoms and the induced surface charges of the continuum solvation model are self‐consistently and variationally determined together with the QM wavefunction. The MM force field methods can be user specified, or a standard force field such as MMFF94, Chemistry at Harvard Molecular Mechanics (CHARMM), Assisted Model Building with Energy Refinement (AMBER), and Optimized Potentials for Liquid Simulations‐All Atom (OPLS‐AA). Analytic gradients for all of these methods are implemented so geometry optimization and molecular dynamics (MD) simulation can be performed. MD free energy perturbation and umbrella sampling methods are also implemented. © 2013 Wiley Periodicals, Inc.     

Monitoring biological membrane-potential changes: a CI QM/MM study

In recent decades, new less-invasive, nonlinear optical methods have been proposed and optimized for monitoring fast physiological processes in biological cells. One of these methods allows the action potential (AP) in cardiomyocytes or neurons to be monitored by means of second-harmonic generation (SHG). We now present the first, to our knowledge, simulations of the dependency of the intensity of the second harmonic (I(SHG)) on variations of the transmembrane potential (TMP) in a cardiomyocyte during an action potential (AP). For this, an amphiphilic potential-sensitive styryl dye molecule with nonlinear optical properties was embedded in a dipalmitoylphosphatidylcholine (DPPC) bilayer, replacing one of the phospholipid molecules. External electrical fields with different strengths were applied across the membrane to simulate the AP of a heart-muscle cell. We used a combined classical/quantum mechanical approach to model the structure and the spectroscopic properties of the embedded chromophore. Two 10 ns molecular dynamics (MD) simulations provided input geometries for semiempirical molecular orbital (QM/MM) single-point configuration interaction (CI) calculations, which were used to calculate the wavelengths and oscillator strengths of electronic transitions in the di-8-ANEPPS dye molecule. The results were then used in a sum-over-states treatment to calculate the second-order hyperpolarizability. The square of the hyperpolarizability scales with the intensity of the second harmonic, which is used to monitor the action potentials of cardiomyocytes experimentally. Thus, we computed changes in the intensity of the second harmonic (DeltaI(SHG)) as function of TMP changes. Our results agree well with experimental measurements.