Effect of protein environment within cytochrome P450cam evaluated using a polarizable-embedding QM/MM method

QM/MM-Based Energy Decomposition Analysis Method for Large Systems

In this work, a QM/MM-based EDA method, called GKS-EDA(QM/MM), is proposed. As an extension of GKS-EDA, this scheme divides the total interaction energy into electrostatic, exchange-repulsion, polarization, and correlation/dispersion terms. GKS-EDA(QM/MM) can be applied to describe the interactions of large-scale systems combined with various QM/MM platforms. By using the examples of a hydrated hydronium ion complex in water solution, the barnase-barstar complex, and MMP-13-pyrimidinetrione in a metalloprotein, the capability of GKS-EDA(QM/MM) for various interactions in large systems is validated.

Multistructural microiteration combined with QM/MM-ONIOM electrostatic embedding

Multistructural microiteration (MSM) is a method to take account of contributions of multiple surrounding structures in a geometrical optimization or reaction path calculation using the quantum mechanics/molecular mechanics (QM/MM) ONIOM method. In this study, we combined MSM with the electrostatic embedding (EE) scheme of the QM/MM-ONIOM method by extending its original formulation for mechanical embedding (ME). MSM-EE takes account of the polarization in the QM region induced by point charges assigned to atoms in the multiple surrounding structures, where the point charges are scaled by the weight factor of each surrounding structure determined through MSM. The performance of MSM-EE was compared with that of the other methods, i.e., ONIOM-ME, ONIOM-EE, and MSM-ME, by applying them to three chemical processes: (1) chorismate-to-prephenate transformation in aqueous solution, (2) the same transformation as (1) in an enzyme, and (3) hydroxylation in p-hydroxybenzoate hydroxylase. These numerical tests of MSM-EE yielded barriers and reaction energies close to experimental values with computational costs comparable to those of the other three methods.

Energy Decomposition Analysis of the Nature of Coordination Bonding at the Heme Iron Center in Cytochrome P450 Inhibition

Drug compounds or their metabolic intermediates (MIs) sometimes inhibit the function of cytochrome P450 enzymes (P450s) by forming a coordination bond with the Fe(III) heme or Fe(II) heme of P450s. Such inhibition is one of the major causes of drug-drug interactions (DDIs), a subject of longstanding academic and practical interest. However, such coordination bonding is not fully understood at the quantum mechanical level, thus hampering rational improvement of the accuracy of DDI-related predictions. In this work, we employed density functional theory (DFT) and the generalized Kohn-Sham energy decomposition analysis (GKS-EDA) scheme to investigate the nature of the coordination bonding formed in the reversible and quasi-irreversible inhibition of P450s. The GKS-EDA results highlighted a previously unrecognized role of the electron correlation effect in P450 inhibition. The correlation effect tends to be larger in Fe(II) complexes of MI-type inhibitors and is particularly prominent for the nitrosoalkane ligand. An additional natural bond orbital (NBO) analysis provided insight into the relative significance of the σ donation and π backdonation effects in various heme-inhibitor complexes.

Harder, better, faster, stronger: Large-scale QM and QM/MM for predictive modeling in enzymes and proteins

Computational prediction of enzyme mechanism and protein function requires accurate physics-based models and suitable sampling. We discuss recent advances in large-scale quantum mechanical (QM) modeling of biochemical systems that have reduced the cost of high-accuracy models. Tradeoffs between sampling and accuracy have motivated modeling with molecular mechanics (MM) in a multiscale QM/MM or iterative approach. Limitations to both conventional density-functional theory and classical MM force fields remain for describing noncovalent interactions in comparison to experiment or wavefunction theory. Because predictions of enzyme action (i.e. electrostatics), free energy barriers, and mechanisms are sensitive to the protocol and embedding method in QM/MM, convergence tests and systematic methods for quantifying QM-level interactions are a needed, active area of development.

Both Configuration and QM Region Size Matter: Zinc Stability in QM/MM Models of DNA Methyltransferase

Quantum-mechanical/molecular-mechanical (QM/MM) methods are essential to the study of metalloproteins, but the relative importance of sampling and degree of QM treatment in achieving quantitative predictions is poorly understood. We study the relative magnitude of configurational and QM-region sensitivity of energetic and electronic properties in a representative Zn²⁺ metal binding site of a DNA methyltransferase. To quantify property variations, we analyze snapshots extracted from 250 ns of molecular dynamics simulation. To understand the degree of QM-region sensitivity, we perform analysis using QM regions ranging from a minimal 49-atom region consisting only of the Zn²⁺ metal and its four coordinating Cys residues up to a 628-atom QM region that includes residues within 12 Å of the metal center. Over the configurations sampled, we observe that illustrative properties (e.g., rigid Zn²⁺ removal energy) exhibit large fluctuations that are well captured with even minimal QM regions. Nevertheless, for both energetic and electronic properties, we observe a slow approach to asymptotic limits with similarly large changes in absolute values that converge only with larger (ca. 300-atom) QM region sizes. For the smaller QM regions, the electronic description of Zn²⁺ binding is incomplete: the metal binds too tightly and is too stabilized by the strong electrostatic potential of MM point charges, and the Zn-S bond covalency is overestimated. Overall, this work suggests that efficient sampling with QM/MM in small QM regions is an effective method to explore the influence of enzyme structure on target properties. At the same time, accurate descriptions of electronic and energetic properties require a larger QM region than the minimal metal-coordinating residues in order to converge treatment of both metal-local bonding and the overall electrostatic environment.

Revealing quantum mechanical effects in enzyme catalysis with large-scale electronic structure simulation

Enzymes have evolved to facilitate challenging reactions at ambient conditions with specificity seldom matched by other catalysts. Computational modeling provides valuable insight into catalytic mechanism, and the large size of enzymes mandates multi-scale, quantum mechanical-molecular mechanical (QM/MM) simulations. Although QM/MM plays an essential role in balancing simulation cost to enable sampling with full QM treatment needed to understand electronic structure in enzyme active sites, the relative importance of these two strategies for understanding enzyme mechanism is not well known. We explore challenges in QM/MM for studying the reactivity and stability of three diverse enzymes: i) Mg2+-dependent catechol O-methyltransferase (COMT), ii) radical enzyme choline trimethylamine lyase (CutC), and iii) DNA methyltransferase (DNMT1), which has structural Zn2+ binding sites. In COMT, strong non-covalent interactions lead to long range coupling of electronic structure properties across the active site, but the more isolated nature of the metallocofactor in DNMT1 leads to faster convergence of some properties. We quantify these effects in COMT by computing covariance matrices of by-residue electronic structure properties during dynamics and along the reaction coordinate. In CutC, we observe spontaneous bond cleavage following initiation events, highlighting the importance of sampling and dynamics. We use electronic structure analysis to quantify the relative importance of CHO and OHO non-covalent interactions in imparting reactivity. These three diverse cases enable us to provide some general recommendations regarding QM/MM simulation of enzymes.

Torsional flexibility of undecorated catechol diether compound as potent NNRTI targeting HIV-1 reverse transcriptase

Conformational adaptation of non-nucleoside reverse transcriptase inhibitor (NNRTI) via torsional flexibility is found to be very significant for targeting human immunodeficiency virus type 1 reverse transcriptase (HIV-1 RT) mutants. Catechol diether derivative including flexible torsions is new potent NNRTI with picomolar activity. Moreover, this derivative also reveals the good solubility, low toxicity and potent inhibition for HIV-1 mutants. In this study, torsional flexibility of an undecorated catechol diether compound in the binding pocket of wild type and mutants (Y181C and K103N/Y181C) HIV-1 RT is investigated by using QM/MM calculations. From the results, the uracil ring is found to exhibit more flexibility in the NNIBP. On the contrary, potential energy surfaces show that high energy is encountered by changing of the corresponding torsion of the cyanovinyl aryl ring indicating the limitation for torsional flexibility. For pointing out the key interaction for the binding, the residual interaction energies are performed by means of QM calculations. Important attractive interactions through hydrogen bonds between the inhibitor and K102, K/N103, V106, and Y188 are observed. The catechol ring is proposed to be modified in order to strengthen interactions with surrounding amino acids. The results may help for the designing of new potent NNRTIs.

Large-scale QM/MM free energy simulations of enzyme catalysis reveal the influence of charge transfer

Hybrid quantum mechanical-molecular mechanical (QM/MM) simulations provide key insights into enzyme structure-function relationships. Numerous studies have demonstrated that large QM regions are needed to systematically converge ground state, zero temperature properties with electrostatic embedding QM/MM. However, it is not well known if ab initio QM/MM free energy simulations have this same dependence, in part due to the hundreds of thousands of energy evaluations required for free energy estimations that in turn limit QM region size. Here, we leverage recent advances in electronic structure efficiency and accuracy to carry out range-separated hybrid density functional theory free energy simulations in a representative methyltransferase. By studying 200 ps of ab initio QM/MM dynamics for each of five QM regions from minimal (64 atoms) to one-sixth of the protein (544 atoms), we identify critical differences between large and small QM region QM/MM in charge transfer between substrates and active site residues as well as in geometric structure and dynamics that coincide with differences in predicted free energy barriers. Distinct geometric and electronic structure features in the largest QM region indicate that important aspects of enzymatic rate enhancement in methyltransferases are identified with large-scale electronic structure.

Revisiting the catalytic mechanism of Mo-Cu carbon monoxide dehydrogenase using QM/MM and DFT calculations

Previous density functional theory (DFT) studies have shown that the release of the produced carbon dioxide (CO2) from an active-site cluster is a thermodynamically or kinetically difficult step in the enzymatic carbon monoxide (CO) oxidation catalyzed by Mo-Cu carbon monoxide dehydrogenase (Mo-Cu CODH). To better understand the effect of the protein environment on this difficult CO2 release step as well as other reaction steps, we applied hybrid quantum mechanics and molecular mechanics (QM/MM) calculations to the Mo-Cu CODH enzyme. The results show that in the first step, the equatorial Mo[double bond, length as m-dash]O group in the active-site cluster attacks the nearby CO molecule bound to the Cu site. Afterward, a stable thiocarbonate intermediate is formed in which the CO2 molecule is embedded and the copper-S(μ-sulfido) bond is broken. A free CO2 molecule, i.e., the final product, is then released from the active-site cluster, not directly from the thiocarbonate intermediate but via a previously formed intermediate that also contains CO2 but retains the Cu-S(μ-sulfido) bond. In contrast to the previous DFT results, the calculated barrier for this process was low in our QM/MM calculations. An additional QM/MM analysis of the barrier height showed that the effect of the protein environment on this barrier lowering is not very large. We found that the reason for the low barrier obtained by QM/MM is that the barrier for CO2 release is already not high at the DFT level. These results allow us to conclude that the CO oxidation reaction passes through the formation of a thiocarbonate intermediate, and that the subsequent CO2 release is kinetically not difficult. Nevertheless, the protein environment has an important role to play in making the latter process thermodynamically favored. No low-barrier pathway for the product release could be obtained for the reaction of n-butylisocyanide, which is consistent with the experimental fact that n-butylisocyanide inhibits Mo-Cu CODH.

Performance of the AMOEBA Water Model in the Vicinity of QM Solutes: A Diagnosis Using Energy Decomposition Analysis

The importance of incorporating solvent polarization effects into the modeling of solvation processes has been well-recognized, and therefore a new generation of hybrid quantum mechanics/molecular mechanics (QM/MM) approaches that accounts for this effect is desirable. We present a fully self-consistent, mutually polarizable QM/MM scheme using the AMOEBA force field, in which the total energy of the system is variationally minimized with respect to both the QM electronic density and the MM induced dipoles. This QM/AMOEBA model is implemented through the Q-Chem/LibEFP code interface and then applied to the evaluation of solute-solvent interaction energies for various systems ranging from the water dimer to neutral and ionic solutes (NH₃, NH₄⁺, CN^-) surrounded by increasing numbers of water molecules (up to 100). In order to analyze the resulting interaction energies, we also utilize an energy decomposition analysis (EDA) scheme which identifies contributions from permanent electrostatics, polarization, and van der Waals (vdW) interaction for the interaction between the QM solute and the solvent molecules described by AMOEBA. This facilitates a component-wise comparison against full QM calculations where the corresponding energy components are obtained via a modified version of the absolutely localized molecular orbitals (ALMO)-EDA. The results show that the present QM/AMOEBA model can yield reasonable solute-solvent interaction energies for neutral and cationic species, while further scrutiny reveals that this accuracy highly relies on the delicate balance between insufficiently favorable permanent electrostatics and softened vdW interaction. For anionic solutes where the charge penetration effect becomes more pronounced, the QM/MM interface turns out to be unbalanced. These results are consistent with and further elucidate our findings in a previous study using a slightly different QM/AMOEBA model ( Dziedzic et al. J. Chem. Phys. 2016 , 145 , 124106 ). The implications of these results for further refinement of this model are also discussed.

How Large Should the QM Region Be in QM/MM Calculations? The Case of Catechol O-Methyltransferase

Hybrid quantum mechanical-molecular mechanical (QM/MM) simulations are widely used in studies of enzymatic catalysis. Until recently, it has been cost prohibitive to determine the asymptotic limit of key energetic and structural properties with respect to increasingly large QM regions. Leveraging recent advances in electronic structure efficiency and accuracy, we investigate catalytic properties in catechol O-methyltransferase, a prototypical methyltransferase critical to human health. Using QM regions ranging in size from reactants-only (64 atoms) to nearly one-third of the entire protein (940 atoms), we show that properties such as the activation energy approach within chemical accuracy of the large-QM asymptotic limits rather slowly, requiring approximately 500-600 atoms if the QM residues are chosen simply by distance from the substrate. This slow approach to asymptotic limit is due to charge transfer from protein residues to the reacting substrates. Our large QM/MM calculations enable identification of charge separation for fragments in the transition state as a key component of enzymatic methyl transfer rate enhancement. We introduce charge shift analysis that reveals the minimum number of protein residues (approximately 11-16 residues or 200-300 atoms for COMT) needed for quantitative agreement with large-QM simulations. The identified residues are not those that would be typically selected using criteria such as chemical intuition or proximity. These results provide a recipe for a more careful determination of QM region sizes in future QM/MM studies of enzymes.

Mechanistic insights into the catalytic reaction of plant allene oxide synthase (pAOS) via QM and QM/MM calculations

QM cluster and QM/MM protein models have been employed to understand aspects of the reaction mechanism of plant allene oxide synthase (pAOS). In this study we have investigated two reaction mechanisms for pAOS. The standard pAOS mechanism was contrasted with an alternative involving an additional active site molecule which has been shown to facilitate proton coupled electron transfer (PCET) in related systems. Firstly, we found that the results from QM/MM protein model are comparable with those from the QM cluster model, presumably due to the large active site used. Furthermore, the results from the QM cluster model show that the Fe(III) and Fe(IV) pathways for the standard mechanism have similar energetic and structural properties, indicating that the reaction mechanism may well proceed via both pathways. However, while the PCET process is facilitated by an additional active site bound water in other related families, in pAOS it is not, suggesting this type of process is not general to all closely related family members.

Applications of density functional theory to iron-containing molecules of bioinorganic interest

The past decades have seen an explosive growth in the application of density functional theory (DFT) methods to molecular systems that are of interest in a variety of scientific fields. Owing to its balanced accuracy and efficiency, DFT plays particularly useful roles in the theoretical investigation of large molecules. Even for biological molecules such as proteins, DFT finds application in the form of, e.g., hybrid quantum mechanics and molecular mechanics (QM/MM), in which DFT may be used as a QM method to describe a higher prioritized region in the system, while a MM force field may be used to describe remaining atoms. Iron-containing molecules are particularly important targets of DFT calculations. From the viewpoint of chemistry, this is mainly because iron is abundant on earth, iron plays powerful (and often enigmatic) roles in enzyme catalysis, and iron thus has the great potential for biomimetic catalysis of chemically difficult transformations. In this paper, we present a brief overview of several recent applications of DFT to iron-containing non-heme synthetic complexes, heme-type cytochrome P450 enzymes, and non-heme iron enzymes, all of which are of particular interest in the field of bioinorganic chemistry. Emphasis will be placed on our own work.