Hybrid Quantum Mechanics-molecular Mechanics Research Articles

Assessment of Embedding Schemes in a Hybrid Machine Learning/Classical Potentials (ML/MM) Approach.

Machine learning (ML) methods have reached high accuracy levels for the prediction of in vacuo molecular properties. However, the simulation of large systems solely through ML methods (such as those based on neural network potentials) is still a challenge. In this context, one of the most promising frameworks for integrating ML schemes in the simulation of complex molecular systems are the so-called ML/MM methods. These multiscale approaches combine ML methods with classical force fields (MM), in the same spirit as the successful hybrid quantum mechanics-molecular mechanics methods (QM/MM). The key issue for such ML/MM methods is an adequate description of the coupling between the region of the system described by ML and the region described at the MM level. In the context of QM/MM schemes, the main ingredient of the interaction is electrostatic, and the state of the art is the so-called electrostatic-embedding. In this study, we analyze the quality of simpler mechanical embedding-based approaches, specifically focusing on their application within a ML/MM framework utilizing atomic partial charges derived in vacuo. Taking as reference electrostatic embedding calculations performed at a QM(DFT)/MM level, we explore different atomic charges schemes, as well as a polarization correction computed using atomic polarizabilites. Our benchmark data set comprises a set of about 80k small organic structures from the ANI-1x and ANI-2x databases, solvated in water. The results suggest that the minimal basis iterative stockholder (MBIS) atomic charges yield the best agreement with the reference coupling energy. Remarkable enhancements are achieved by including a simple polarization correction.

Binding mechanism of disulfide species to ferric hemeproteins: The case of metmyoglobin

The interactions of the heme iron of hemeproteins with sulfide and disulfide compounds are of potential interest as physiological signaling processes. While the interaction with hydrogen sulfide has been described computationally and experimentally, the reaction with disulfide, and specifically the molecular mechanism for ligand binding has not been studied in detail. In this work, we study the association process for disulfane and its conjugate base disulfanide at different pH conditions. Additionally, by means of advanced sampling techniques based on multiple steered molecular dynamics, we provide free energy profiles for ligand migration for both acid/base species, showing a similar behavior to the previously reported for the related H2S/HS¯ pair. Finally, we studied the ligand interchange reaction (H2O/H2S, HS¯ and H2O/HSSH, HSS¯) by means of hybrid quantum mechanics-molecular mechanics calculations. We show that the anionic species are able to displace more efficiently the H2O bound to the iron, and that the H-bond network in the distal cavity can help the neutral species to perform the reaction. Altogether, we provide a molecular explanation for the experimental information and show that the global association process depends on a fine balance between the migration towards the active site and the ligand interchange reaction.

Multiscale Modeling Unravels the Influence of Biomembranes on the Photochemical Properties of Embedded Anti-Oxidative Polyphenolic and Phenanthroline Chelating Dyes.

The embedding of caffeate methyl ester, the flavonoids luteolin and quercetin, and the o-phenanthroline and neocuproine in a liquid disordered lipid bilayer has been studied through extensive atomistic calculations. The location and the orientation of these bio-active antioxidants are explained and analyzed. While the two phenanthrolines strongly associate with the lipid tail region, the other three compounds are rather found among the head groups. The simulations showcase conformational changes of the flavonoids. Through the use of a hybrid quantum mechanics-molecular mechanics scheme and supported by a profound benchmarking of the electronic excited-state method for these compounds, the influence of the anisotropic environment on the compounds' optical properties is analyzed. Influences of surrounding water molecules and of the polar parts of the lipids on the transition dipole moments and excited-state dipole moments are weighted with respect to a change in conformation. The current study highlights the importance of the mapping of molecular interactions in model membranes and pinpoints properties, which can be biomedically used to discriminate and detect different lipid environments.

A Benchmark Study of Quantum Mechanics and Quantum Mechanics-Molecular Mechanics Methods for Carbocation Chemistry

Carbocations play key roles in classical organic reactions and have also been implicated in several enzyme families. A hallmark of carbocation chemistry is multitudes of competing reaction pathways, and to be able to distinguish between pathways with quantum chemical calculations, it is necessary to approach chemical accuracy for relative energies between carbocations. Here, we present an extensive study of the performance of selected density functional theory (DFT) methods in describing the thermochemistry and kinetics of carbocations and their corresponding neutral alkenes both in the gas-phase and within a hybrid quantum mechanics-molecular mechanics (QM/MM) framework. The density functionals are benchmarked against accurate ab initio methods such as CBS-QB3 and DLPNO-CCSD(T). Based on the findings in the gas-phase calculations of carbocations and alkenes, the best functionals are chosen and tested further for non-covalent interactions in model systems using QM and QM/MM methods. We compute the interaction energies between a model carbocation/alkane and model π, dipole, and hydrophobic systems using DFT and QM(DFT)/MM and compare with DLPNO-CCSD(T). These latter model systems are representative of side chains of amino acids such as phenylalanine/tyrosine, tryptophan, asparagine/glutamine, serine/threonine, methionine, and other hydrophobic groups. The Lennard-Jones parameters of the QM atoms in QM(DFT)/MM calculations are modified to obtain an optimal fit with the QM energies. Finally, a selected carbocation reaction is studied in the gas phase and in implicit chloroform solvent using QM and in explicit chloroform solvent using QM/MM and umbrella sampling simulations. This study highlights the highest accuracy possible with selected density functionals and QM/MM methods but also some limitations in using QM/MM methods for carbocation systems.

Mechanism of Nucleotidyltransfer Reaction and Role of Mg2+ Ion in Sugar Nucleotidyltransferases

Sugar nucleotidyl transferases (SNTs) catalyze nucleotidyltransfer reactions to form sugar-nucleotides and pyrophosphate in the presence of two Mg2+ ions (Mg2+A and Mg2+B). We unveil the mechanism and free energetics of nucleotidyl transfer reaction in an SNT called GlmU through hybrid quantum mechanics-molecular mechanics molecular dynamics simulations and free energy calculations. The study identifies the roles of the active site residues and the Mg2+ ions in catalyzing the reaction. Of great significance, we are able to compare the free energy barrier for the reaction with that for the Mg2+-assisted release of the product (i.e., pyrophosphate) into the solution, shedding light on the general mechanistic and kinetic aspects of catalysis by SNTs.

EnzyDock: Protein–Ligand Docking of Multiple Reactive States along a Reaction Coordinate in Enzymes

Enzymes play a pivotal role in all biological systems. These biomachines are the most effective catalysts known, dramatically enhancing the rate of reactions by more than 10 orders of magnitude relative to the uncatalyzed reactions in solution. Predicting the correct, mechanistically appropriate binding modes for substrate and product, as well as all reaction intermediates and transition states, along a reaction pathway is immensely challenging and remains an unsolved problem. In the present work, we developed an effective methodology for identifying probable binding modes of multiple ligand states along a reaction coordinate in an enzyme active site. The program is called EnzyDock and is a CHARMM-based multistate consensus docking program that includes a series of protocols to predict the chemically relevant orientation of substrate, reaction intermediates, transition states, product, and inhibitors. EnzyDock is based on simulated annealing molecular dynamics and Monte Carlo sampling and allows ligand, as well as enzyme side-chain and backbone flexibility. The program can employ many user-defined constraints and restraints and classical force field potentials, as well as a range of hybrid quantum mechanics-molecular mechanics potentials. Herein, we apply EnzyDock to several different kinds of problems. First, we study two terpene synthase reactions, namely bornyl diphosphate synthase and the bacterial diterpene synthase CotB2. Second, we use EnzyDock to predict reaction coordinate states in a pair of Diels-Alder reactions in the enzymes spirotetronate AbyU and LepI. Third, we study a couple of racemases: the cofactor-dependent serine racemase and the cofactor independent proline racemase. Finally, we study several cases of covalent docking involving the Michael addition reaction. For all systems we predict binding modes that are consistent with available experimental observations, as well as with theoretical modeling studies from the literature. EnzyDock provides a platform for generating mechanistic insight into enzyme reactions, useful and reliable starting points for in-depth multiscale modeling projects, and rational design of noncovalent and covalent enzyme inhibitors.

Slow-Starter Enzymes: Role of Active-Site Architecture in the Catalytic Control of the Biosynthesis of Taxadiene by Taxadiene Synthase

Taxadiene synthase (TXS) catalyzes the formation of natural product taxa-4(5),11(12)-diene (henceforth taxadiene). Taxadiene is the precursor in the formation of Taxol, which is an important natural anticancer agent. In the current study, we present a detailed mechanistic view of the biosynthesis of taxadiene by TXS using a hybrid quantum mechanics-molecular mechanics potential in conjunction with free energy simulation methods. The obtained free-energy landscape displays initial endergonic steps followed by a stepwise downhill profile, which is an emerging free-energy fingerprint for type I terpene synthases. We identify an active-site Trp residue (W753) as a key feature of the TXS active-site architecture and propose that this residue stabilized intermediate cations via π-cation interactions. To validate our proposed active TXS model, we examine a previously reported W753H mutation, which leads to the exclusive formation of side product cembrene A. The simulations of the W753H mutant show that, in the mutant structure, the His side chain is in the perfect position to deprotonate the cembrenyl cation en route to cembrene formation and that this abortive deprotonation is an energetically facile process. On the basis of the current model, we propose that an analogous mutation of Y841 to His could possibly lead to verticillane. The current simulations stress the importance of the precise positioning of key active-site residues in stabilizing intermediate carbocations. In view of the great pharmaceutical importance of taxadiene, a detailed understanding of the TXS mechanism can provide important clues toward a synthetic strategy for Taxol manufacturing.

The role of the Met20 loop in the hydride transfer in Escherichia coli dihydrofolate reductase

A key question concerning the catalytic cycle of Escherichia coli dihydrofolate reductase (ecDHFR) is whether the Met20 loop is dynamically coupled to the chemical step during catalysis. A more basic, yet unanswered question is whether the Met20 loop adopts a closed conformation during the chemical hydride transfer step. To examine the most likely conformation of the Met20 loop during the chemical step, we studied the hydride transfer in wild type (WT) ecDHFR using hybrid quantum mechanics-molecular mechanics free energy simulations with the Met20 loop in a closed and disordered conformation. Additionally, we investigated three mutant forms (I14X; X = Val, Ala, Gly) of the enzyme that have increased active site flexibility and donor-acceptor distance dynamics in closed and disordered Met20 loop states. We found that the conformation of the Met20 loop has a dramatic effect on the ordering of active site hydration, although the Met20 loop conformation only has a moderate effect on the hydride transfer rate and donor-acceptor distance dynamics. Finally, we evaluated the pKa of the substrate N5 position in closed and disordered Met20 loop states and found a strong correlation between N5 basicity and the conformation of the Met20 loop.

Improved Sugar Puckering Profiles for Nicotinamide Ribonucleoside for Hybrid QM/MM Simulations.

The coenzyme nicotinamide adenine dinucleotide (NAD+) and its reduced form (NADH) play ubiquitous roles as oxidizing and reducing agents in nature. The binding, and possibly the chemical redox step, of NAD+/NADH may be influenced by the cofactor conformational distribution and, in particular, by the ribose puckering of its nicotinamide-ribonucleoside (NR) moiety. In many hybrid quantum mechanics-molecular mechanics (QM/MM) studies of NAD+/NADH dependent enzymes, the QM region is treated by semiempirical (SE) methods. Recent work suggests that SE methods do not adequately describe the ring puckering in sugar molecules. In the present work we adopt an efficient and practical strategy to correct for this deficiency for NAD+/NADH. We have implemented a cost-effective correction to a SE Hamiltonian by adding a correction potential, which is defined as the difference between an accurate benchmark density functional theory (DFT) potential energy surface (PES) and the SE PES. In practice, this is implemented via a B-spline interpolation scheme for the grid-based potential energy difference surface. We find that the puckering population distributions obtained from free energy QM(SE)/MM simulations are in good agreement with DFT and in fair accord with experimental results. The corrected PES should facilitate a more accurate description of the ribose puckering in the NAD+/NADH cofactor in simulations of biological systems.

First principles model calculations of the biosynthetic pathway in selinadiene synthase

Terpenes comprise the largest class of natural products currently known. These ubiquitous molecules are synthesized by terpene synthases via complex carbocationic reactions, incorporating highly reactive intermediates. In the current study, we present a mechanistic investigation of the biosynthetic pathway for the formation of selina-4(15),7(11)-diene. We employ density functional theory to study a model carbocation system in the gas-phase, and delineate the energetic feasibility of a plausible reaction path. Our results suggests that during formation of selina-4(15),7(11)-diene, the substrate is likely folded in a conformation conducive to sequential cyclizations. We propose that a required proton transfer cannot occur intramolecularly in the gas-phase due to a high free energy barrier, and that enzyme assistance is essential for this step. Hybrid quantum mechanics-molecular mechanics docking studies suggest that enzyme intervention could be realized through electrostatic guidance.

Long-range electrostatic corrections in multipolar/polarizable QM/MM simulations.

Taking long-range electrostatic effects into account in classical and hybrid quantum mechanics-molecular mechanics (QM/MM) simulations is necessary for an accurate description of the system under study. We have recently developed a method, termed long-range electrostatic corrections (LREC), for monopolar QM/MM calculations. Here, we present an extension of LREC for multipolar/polarizable QM/MM simulations within the LICHEM software package. Reaction barriers and QM-MM interaction energies converge with a LREC cutoff between 20 and 25 Å, in agreement with our previous results. Additionally, the LREC approach for the QM-MM interactions can be smoothly combined with standard shifting or Ewald summation methods in the MM calculations. We recommend the use of QM(LREC)/MM(PME), where the QM region is treated with LREC and the MM region is treated with particle mesh Ewald (PME) summation. This combination is an excellent compromise between simplicity, speed, and accuracy for large QM/MM simulations.

One- and two-photon absorption of a spiropyran-merocyanine system: experimental and theoretical studies.

We report on the nonlinear optical properties measurements and quantum-chemical calculations of a well-known photochromic system consisting of spiropyran and the merocyanine photoproduct. The study of nonlinear absorption and refraction properties of the molecules dissolved in chloroform were performed with the Z-scan technique, using femtosecond pulses in a wide range of wavelengths. Maxima in the two-photon absorption spectrum at 700 and 1050 nm were found for the merocyanine form, and the corresponding two-photon absorption cross section is 80 GM and 20 GM, respectively. The latter feature does not vanish completely in the nonlinear spectrum of the spiropyran form, possibly because of the existence of some photoconversion caused by the laser beam during the measurements. A nonlinear absorption peak at 900 nm is found in the spiropyran form with an effective cross section of about 20 GM; it is likely due to three-photon absorption or to absorption by some intermediate species. The experimental data are supported by calculations performed with the use of a hybrid quantum mechanics-molecular mechanics approach.

Simulations of remote mutants of dihydrofolate reductase reveal the nature of a network of residues coupled to hydride transfer.

Recent experimental and theoretical studies have proposed that enzymes involve networks of coupled residues throughout the protein that participate in motions accompanying chemical barrier crossing. Here, we have examined portions of a proposed network in dihydrofolate reductase (DHFR) using quantum mechanics/molecular mechanics simulations. The simulations use a hybrid quantum mechanics-molecular mechanics approach with a recently developed semiempirical AM1-SRP Hamiltonian that provides accurate results for this reaction. The simulations reproduce experimentally determined catalytic rates for the wild type and distant mutants of E. coli DHFR, underscoring the accuracy of the simulation protocol. Additionally, the simulations provide detailed insight into how residues remote from the active site affect the catalyzed chemistry, through changes in the thermally averaged properties along the reaction coordinate. The mutations do not greatly affect the structure of the transition state near the bond activation, but we observe differences somewhat removed from the point of C-H cleavage that affect the rate. The mutations have global effects on the thermally averaged structure that propagate throughout the enzyme and the current simulations highlight several interactions that appear to be particularly important.

Electric dipole moments of the fluorescent probes Prodan and Laurdan: experimental and theoretical evaluations.

Several experimental and theoretical approaches can be used for a comprehensive understanding of solvent effects on the electronic structure of solutes. In this review, we revisit the influence of solvents on the electronic structure of the fluorescent probes Prodan and Laurdan, focusing on their electric dipole moments. These biologically used probes were synthesized to be sensitive to the environment polarity. However, their solvent-dependent electronic structures are still a matter of discussion in the literature. The absorption and emission spectra of Prodan and Laurdan in different solvents indicate that the two probes have very similar electronic structures in both the ground and excited states. Theoretical calculations confirm that their electronic ground states are very much alike. In this review, we discuss the electric dipole moments of the ground and excited states calculated using the widely applied Lippert-Mataga equation, using both spherical and spheroid prolate cavities for the solute. The dimensions of the cavity were found to be crucial for the calculated dipole moments. These values are compared to those obtained by quantum mechanics calculations, considering Prodan in vacuum, in a polarizable continuum solvent, and using a hybrid quantum mechanics-molecular mechanics methodology. Based on the theoretical approaches it is evident that the Prodan dipole moment can change even in the absence of solute-solvent-specific interactions, which is not taken into consideration with the experimental Lippert-Mataga method. Moreover, in water, for electric dipole moment calculations, it is fundamental to consider hydrogen-bonded molecules.

Quantum and Classical Simulations of Orotidine Monophosphate Decarboxylase: Support for a Direct Decarboxylation Mechanism

Orotidine 5'-monophosphate (OMP) decarboxylase (ODCase) catalyzes the decarboxylation of OMP to uridine 5'-monophosphate (UMP). Numerous studies of this reaction have suggested a plethora of mechanisms including covalent addition, ylide or carbene formation, and concerted or stepwise protonation. Recent experiments and simulations present strong evidence for a direct decarboxylation mechanism, although direct comparison between experiment and theory is still lacking. In the current work we present hybrid quantum mechanics-molecular mechanics simulations that address the detailed decarboxylation mechanisms for OMP and 5-fluoro-OMP by ODCase. Multidimensional potentials of mean force are computed as functions of structural progress coordinates for the Methanobacterium thermoautotrophicum ODCase reaction: the decarboxylation reaction coordinate, an orbital rehybridization coordinate, and the proton transfer coordinate between Lys72 and the substrate. The computed free energy profiles are in accord with the available experimental data. To facilitate further direct comparison with experiment, we compute the kinetic isotope effects (KIEs) for the enzyme-catalyzed reactions using a mass-perturbation-based path-integral method. The computed KIE provide further support for a direct decarboxylation mechanism. In agreement with experiment, the data suggest a role for Lys72 in stabilizing the transition state in the catalysis of OMP and, to a somewhat lesser extent, in 5-fluoro-OMP.

Structure of Carboxyl-Acid-Terminated Self-Assembled Monolayers from Molecular Dynamics Simulations and Hybrid Quantum Mechanics–Molecular Mechanics Vibrational Normal-Mode Analysis

Self-assembled monolayers (SAMs) are excellent models for studying interfacial reactions and various properties of thin films. Carboxylic acid-terminated SAMs (CATSAMs) are of special interest because their surface properties are highly pH-dependent. To elucidate the complicated pH-dependent structural properties of CATSAMs, we combine force-field molecular dynamics simulations with hybrid quantum mechanics-molecular mechanics (QM/MM) calculations of vibrational frequencies. These studies show that at low pH the carboxylic acids behave similarly to bulk carboxylic acids, hydrogen bonding mainly to the aqueous phase. With increasing pH, intralayer hydrogen bonds and monolayer disorder increase. The computed QM/MM vibrational data show that the observed range of carbonyl frequencies is due to inherent condensed phase properties of the CATSAMs such as fluctuating hydrogen bond environment and steric interactions. The results reported herein underscore the role H-bonding and steric congestion play in pH-depen...

Solvatochromic shift of phenol blue in water from a combined Car–Parrinello molecular dynamics hybrid quantum mechanics-molecular mechanics and ZINDO approach

The present work addresses the solvatochromic shift of phenol blue (PB) dye. For this purpose the results of Car-Parrinello molecular dynamics (CPMD) simulations for PB in gas phase are compared with results obtained for PB in water from CPMD hybrid quantum mechanics-molecular mechanics (CPMD-QM/MM) calculations. The absorption spectra were obtained using the intermediate neglect of differential overlap/spectroscopic-configuration interaction (INDO/CIS) method and were calculated for a multitude of configurations of the trajectory. The calculated lambda(max) for PB in gas phase was found to be about 535 nm, which is considerably lower than the lambda(max) reported for PB in nonpolar solvents. Different solvation shells for PB in water have been defined based on the solute-all-atoms and solvent center of mass radial distribution function (g(r(X-O))). The electronic excitation energies for PB computed in the presence of solvent molecules in an increasing number of solvation shells were calculated in a systematic way to evaluate their contributions to the solvatochrmic shift. The inclusion of solvent molecules in the hydration shell yields a lambda(max) of 640 nm, which contributes to almost 78% of the solvatochromic shift. The inclusion of solvent molecules up to 10 A in the g(r(X-O)) rdf yields a lambda(max) of 670 nm which is in good agreement with the experimentally reported value of 654-684 nm. Overall, the present study suggests that the combined CPMD-QM/MM and INDO-CIS approach can be used successfully to model solvatochromic shifts of organic dye molecules.

What Governs the Charge Transfer in DNA? The Role of DNA Conformation and Environment

Charge transfer in DNA has received much attention in the last few years due to its role in oxidative damage and repair in DNA and also due to possible applications of DNA in nanoelectronics. Despite intense experimental and theoretical efforts, the mechanism underlying long-range hole transport is still unresolved. This is in particular due to the sensitive dependence of charge transfer on the complex structure and dynamics of DNA and the interaction with the solvent, which could not be addressed adequately in the modeling approaches up to now. In this work, we study the factors governing hole transfer in detail, using a DFT-based fragment-orbital method, which allows to compute the charge transfer parameters along multinanosecond molecular dynamics simulations. Environmental effects are captured using a hybrid quantum mechanics-molecular mechanics (QM/MM) coupling scheme. This methodology allows to analyze several factors responsible for charge transfer in DNA in detail. The fluctuation of counterions, strongly counterbalanced by the surrounding water, leads to large oscillations of onsite energies, which govern the energetics of hole propagation along the DNA strand. In contrast, the electronic couplings depend only on DNA conformation and are not affected by the solvent. In particular, the onsite energies are strongly correlated between neighboring nucleobases, indicating that a conformational-gating type of mechanism may be induced by the collective environmental degrees of freedom.

Modeling heme proteins using atomistic simulations

Heme proteins are found in all living organisms, and perform a wide variety of tasks ranging from electron transport, to the oxidation of organic compounds, to the sensing and transport of small molecules. In this work we review the application of classical and quantum-mechanical atomistic simulation tools to the investigation of several relevant issues in heme proteins chemistry: (i) conformational analysis, ligand migration, and solvation effects studied using classical molecular dynamics simulations; (ii) electronic structure and spin state energetics of the active sites explored using quantum-mechanics (QM) methods; (iii) the interaction of heme proteins with small ligands studied through hybrid quantum mechanics-molecular mechanics (QM-MM) techniques; (iv) and finally chemical reactivity and catalysis tackled by a combination of quantum and classical tools.