LICHEM 1.1: Recent Improvements and New Capabilities

A Polarizable QM/MM Model That Combines the State-Averaged CASSCF and AMOEBA Force Field for Photoreactions in Proteins

This study presents the polarizable quantum mechanics/molecular mechanics (QM/MM) embedding of the state-averaged complete active space self-consistent field (SA-CASSCF) in the atomic multipole optimized energetics for biomolecular applications (AMOEBA) force field for the purpose of studying photoreactions in protein environments. We describe two extensions of our previous work that combine SA-CASSCF with AMOEBA water models, allowing it to be generalized to AMOEBA models for proteins and other macromolecules. First, we discuss how our QM/MM model accounts for the discrepancy between the direct and polarization electric fields that arises in the AMOEBA description of intramolecular polarization. A second improvement is the incorporation of link atom schemes to treat instances in which the QM/MM boundary goes through covalent bonds. A single-link atom scheme and double-link atom scheme are considered in this work, and we will discuss how electrostatic interaction, van der Waals interaction, and various kinds of valence terms are treated across the boundary. To test the accuracy of the link atom scheme, we will compare QM/MM with full QM calculations and study how the errors in ground state properties, excited state properties, and excitation energies change when tuning the parameters in the link atom scheme. We will also test the new SA-CASSCF/AMOEBA method on an elementary reaction step in NanoLuc, an artificial bioluminescence luciferase. We will show how the reaction mechanism is different when calculated in the gas phase, in polarizable continuum medium (PCM), versus in protein AMOEBA models.

PyDFT-QMMM: A modular, extensible software framework for DFT-based QM/MM molecular dynamics

PyDFT-QMMM is a Python-based package for performing hybrid quantum mechanics/molecular mechanics (QM/MM) simulations at the density functional level of theory. The program is designed to treat short-range and long-range interactions through user-specified combinations of electrostatic and mechanical embedding procedures within periodic simulation domains, providing necessary interfaces to external quantum chemistry and molecular dynamics software. To enable direct embedding of long-range electrostatics in periodic systems, we have derived and implemented force terms for our previously described QM/MM/PME approach [Pederson and McDaniel, J. Chem. Phys. 156, 174105 (2022)]. Communication with external software packages Psi4 and OpenMM is facilitated through Python application programming interfaces (APIs). The core library contains basic utilities for running QM/MM molecular dynamics simulations, and plug-in entry-points are provided for users to implement custom energy/force calculation and integration routines, within an extensible architecture. The user interacts with PyDFT-QMMM primarily through its Python API, allowing for complex workflow development with Python scripting, for example, interfacing with PLUMED for free energy simulations. We provide benchmarks of forces and energy conservation for the QM/MM/PME and alternative QM/MM electrostatic embedding approaches. We further demonstrate a simple example use case for water solute in a water solvent system, for which radial distribution functions are computed from 100 ps QM/MM simulations; in this example, we highlight how the solvation structure is sensitive to different basis-set choices due to under- or over-polarization of the QM water molecule's electron density.

MiMiC: A high-performance framework for multiscale molecular dynamics simulations

MiMiC is a framework for performing multiscale simulations in which loosely coupled external programs describe individual subsystems at different resolutions and levels of theory. To make it highly efficient and flexible, we adopt an interoperable approach based on a multiple-program multiple-data (MPMD) paradigm, serving as an intermediary responsible for fast data exchange and interactions between the subsystems. The main goal of MiMiC is to avoid interfering with the underlying parallelization of the external programs, including the operability on hybrid architectures (e.g., CPU/GPU), and keep their setup and execution as close as possible to the original. At the moment, MiMiC offers an efficient implementation of electrostatic embedding quantum mechanics/molecular mechanics (QM/MM) that has demonstrated unprecedented parallel scaling in simulations of large biomolecules using CPMD and GROMACS as QM and MM engines, respectively. However, as it is designed for high flexibility with general multiscale models in mind, it can be straightforwardly extended beyond QM/MM. In this article, we illustrate the software design and the features of the framework, which make it a compelling choice for multiscale simulations in the upcoming era of exascale high-performance computing.

Multiscale biomolecular simulations in the exascale era

The complexity of biological systems and processes, spanning molecular to macroscopic scales, necessitates the use of multiscale simulations to get a comprehensive understanding. Quantum mechanics/molecular mechanics (QM/MM) molecular dynamics (MD) simulations are crucial for capturing processes beyond the reach of classical MD simulations. The advent of exascale computing offers unprecedented opportunities for scientific exploration, not least within life sciences, where simulations are essential to unravel intricate molecular mechanisms underlying biological processes. However, leveraging the immense computational power of exascale computing requires innovative algorithms and software designs. In this context, we discuss the current status and future prospects of multiscale biomolecular simulations on exascale supercomputers with a focus on QM/MM MD. We highlight our own efforts in developing a versatile and high-performance multiscale simulation framework with the aim of efficient utilization of state-of-the-art supercomputers. We showcase its application in uncovering complex biological mechanisms and its potential for leveraging exascale computing.

Effects of the Y432S Cancer-Associated Variant on the Reaction Mechanism of Human DNA Polymerase κ

Human DNA polymerases are vital for genetic information management. Their function involves catalyzing the synthesis of DNA strands with unparalleled accuracy, which ensures the fidelity and stability of the human genomic blueprint. Several disease-associated mutations and their functional impact on DNA polymerases have been reported. One particular polymerase, human DNA polymerase kappa (Pol κ), has been reported to be susceptible to several cancer-associated mutations. The Y432S mutation in Pol κ, associated with various cancers, is of interest due to its impact on polymerization activity and markedly reduced thermal stability. Here, we have used computational simulations to investigate the functional consequences of the Y432S using classical molecular dynamics (MD) and coupled quantum mechanics/molecular mechanics (QM/MM) methods. Our findings suggest that Y432S induces structural alterations in domains responsible for nucleotide addition and ternary complex stabilization while retaining structural features consistent with possible catalysis in the active site. Calculations of the minimum energy path associated with the reaction mechanism of the wild type (WT) and Y432S Pol κ indicate that, while both enzymes are catalytically competent (in terms of energetics and the active site's geometries), the cancer mutation results in an endoergic reaction and an increase in the catalytic barrier. Interactions with a third magnesium ion and environmental effects on nonbonded interactions, particularly involving key residues, contribute to the kinetic and thermodynamic distinctions between the WT and mutant during the catalytic reaction. The energetics and electronic findings suggest that active site residues favor the catalytic reaction with dCTP3- over dCTP4-.

Impact of an Ionic Liquid Solution on Horseradish Peroxidase Activity

Horseradish peroxidase (HRP) is an enzyme that oxidizes pollutants from wastewater. A previous report indicated that peroxidases can have an enhancement in initial enzymatic activity in an aqueous solution of 0.26 M 1-ethyl-3-methylimidazolium ethyl sulfate ([EMIm][EtSO4]) at neutral pH. However, the atomistic details remain elusive. In the enzymatic landscape of HRP, compound II (Cpd II) plays a key role and involves a histidine (H42) residue. Cpd II exists as oxoferryl (2a) or hydroxoferryl (2b(FeIV)) forms, where 2a is the predominantly observed form in experimental studies. Intriguingly, the ferric 2b(FeIII) form seen in synthetic complexes has not been observed in HRP. Here, we have investigated the structure and dynamics of HRP in pure water and aqueous [EMIm][EtSO4] (0.26 M), as well as the reaction mechanism of 2a to 2b conversion using polarizable molecular dynamics (MD) simulations and quantum mechanics/molecular mechanics (QM/MM) calculations. When HRP is solvated in aq [EMIm][EtSO4], the catalytic water displaces, and H42 directly orients over the ferryl moiety, allowing a direct proton transfer (PT) with a significant energy barrier reduction. Conversely, in neat water, the reaction of 2a to 2b follows the previously reported mechanism. We further investigated the deprotonated form of H42. Analysis of the electric fields at the active site indicates that the aq [EMIm][EtSO4] medium facilitates the reaction by providing a more favorable environment compared with the system solvated in neat water. Overall, the atomic level supports the previous experimental observations and underscores the importance of favorable electric fields in the active site to promote catalysis.

Seamless integration of GEM, a density based-force field, for QM/MM simulations via LICHEM, Psi4, and Tinker-HP

Hybrid quantum mechanics/molecular mechanics (QM/MM) simulations have become an essential tool in computational chemistry, particularly for analyzing complex biological and condensed phase systems. Building on this foundation, our work presents a novel implementation of the Gaussian Electrostatic Model (GEM), a polarizable density-based force field, within the QM/MM framework. This advancement provides seamless integration, enabling efficient and optimized QM/GEM calculations in a single step using the LICHEM Code. We have successfully applied our implementation to water dimers and hexamers, demonstrating the ability to handle water systems with varying numbers of water molecules. Moreover, we have extended the application to describe the double proton transfer of the aspartic acid dimer in a box of water, which highlights the method's proficiency in investigating heterogeneous systems. Our implementation offers the flexibility to perform on-the-fly density fitting or to utilize pre-fitted coefficients to estimate exchange and Coulomb contributions. This flexibility enhances efficiency and accuracy in modeling molecular interactions, especially in systems where polarization effects are significant.

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.

Relative cooperativity in neutral and charged molecular clusters using QM/MM calculations

QM/MM methods have been used to study electronic structure properties and chemical reactivity in complex molecular systems where direct electronic structure calculations are not feasible. In our previous work, we showed that non-polarizable force fields, by design, describe intermolecular interactions through pairwise interactions, overlooking many-body interactions involving three or more particles. In contrast, polarizable force fields account partially for many-body effects through polarization, but still handle van der Waals and permanent electrostatic interactions pairwise. We showed that despite those limitations, polarizable and non-polarizable force fields can reproduce relative cooperativity achieved using density functional theory due to error compensation mechanisms. In this contribution, we assess the performance of QM/MM methods in reproducing these phenomena. Our study highlights the significance of the QM region size and force field choice in QM/MM calculations, emphasizing the importance of parameter validation to obtain accurate interaction energy predictions.

Leveraging QM/MM and Molecular Dynamics Simulations to Decipher the Reaction Mechanism of the Cas9 HNH Domain to Investigate Off-Target Effects

The clustered regularly interspaced short palindromic repeats (CRISPR) technology is an RNA-guided targeted genome-editing tool using Cas family proteins. Two magnesium-dependent nuclease domains of the Cas9 enzyme, termed HNH and RuvC, are responsible for cleaving the target DNA (t-DNA) and nontarget DNA strands, respectively. The HNH domain is believed to determine the DNA cleavage activity of both endonuclease domains and is sensitive to complementary RNA-DNA base pairing. However, the underlying molecular mechanisms of CRISPR-Cas9, by which it rebukes or accepts mismatches, are poorly understood. Thus, investigation of the structure and dynamics of the catalytic state of Cas9 with either matched or mismatched t-DNA can provide insights into improving its specificity by reducing off-target cleavages. Here, we focus on a recently discovered catalytic-active form of the Streptococcus pyogenes Cas9 (SpCas9) and employ classical molecular dynamics and coupled quantum mechanics/molecular mechanics simulations to study two possible mechanisms of t-DNA cleavage reaction catalyzed by the HNH domain. Moreover, by designing a mismatched t-DNA structure called MM5 (C to G at the fifth position from the protospacer adjacent motif region), the impact of single-guide RNA (sgRNA) and t-DNA complementarity on the catalysis process was investigated. Based on these simulations, our calculated binding affinities, minimum energy paths, and analysis of catalytically important residues provide atomic-level details of the differences between matched and mismatched cleavage reactions. In addition, several residues exhibit significant differences in their catalytic roles for the two studied systems, including K253, K263, R820, K896, and K913.

Investigation of the Inhibition Mechanism of Xanthine Oxidoreductase by Oxipurinol: A Computational Study

Xanthine oxidoreductase (XOR) is an enzyme found in various organisms. It converts hypoxanthine to xanthine and urate, which are crucial steps in purine elimination in humans. Elevated uric acid levels can lead to conditions like gout and hyperuricemia. Therefore, there is significant interest in developing drugs that target XOR for treating these conditions and other diseases. Oxipurinol, an analogue of xanthine, is a well-known inhibitor of XOR. Crystallographic studies have revealed that oxipurinol directly binds to the molybdenum cofactor (MoCo) in XOR. However, the precise details of the inhibition mechanism are still unclear, which would be valuable for designing more effective drugs with similar inhibitory functions. In this study, molecular dynamics and quantum mechanics/molecular mechanics calculations are employed to investigate the inhibition mechanism of XOR by oxipurinol. The study examines the structural and dynamic effects of oxipurinol on the pre-catalytic structure of the metabolite-bound system. Our results provide insights on the reaction mechanism catalyzed by the MoCo center in the active site, which aligns well with experimental findings. Furthermore, the results provide insights into the residues surrounding the active site and propose an alternative mechanism for developing alternative covalent inhibitors.

Computational Characterization of the Inhibition Mechanism of Xanthine Oxidoreductase by Topiroxostat

Xanthine oxidase (XO) is a member of the molybdopterin-containing enzyme family. It interconverts xanthine to uric acid as the last step of purine catabolism in the human body. The high uric acid concentration in the blood directly leads to human diseases like gout and hyperuricemia. Therefore, drugs that inhibit the biosynthesis of uric acid by human XO have been clinically used for many years to decrease the concentration of uric acid in the blood. In this study, the inhibition mechanism of XO and a new promising drug, topiroxostat (code: FYX-051), is investigated by employing molecular dynamics (MD) and quantum mechanics/molecular mechanics (QM/MM) calculations. This drug has been reported to act as both a noncovalent and covalent inhibitor and undergoes a stepwise inhibition by all its hydroxylated metabolites, which include 2-hydroxy-FYX-051, dihydroxy-FYX-051, and trihydroxy-FYX-051. However, the detailed mechanism of inhibition of each metabolite remains elusive and can be useful for designing more effective drugs with similar inhibition functions. Hence, herein we present the computational investigation of the structural and dynamical effects of FYX-051 and the calculated reaction mechanism for all of the oxidation steps catalyzed by the molybdopterin center in the active site. Calculated results for the proposed reaction mechanisms for each metabolite's inhibition reaction in the enzyme's active site, binding affinities, and the noncovalent interactions with the surrounding amino acid residues are consistent with previously reported experimental findings. Analysis of the noncovalent interactions via energy decomposition analysis (EDA) and noncovalent interaction (NCI) techniques suggests that residues L648, K771, E802, R839, L873, R880, R912, F914, F1009, L1014, and A1079 can be used as key interacting residues for further hybrid-type inhibitor development.

The application of QM/MM simulations in heterogeneous catalysis

The QM/MM simulation method is provenly efficient for the simulation of biological systems, where an interplay of extensive environment and delicate local interactions drives a process of interest through a funnel on a complex energy landscape. Recent advances in quantum chemistry and force-field methods present opportunities for the adoption of QM/MM to simulate heterogeneous catalytic processes, and their related systems, where similar intricacies exist on the energy landscape. Herein, the fundamental theoretical considerations for performing QM/MM simulations, and the practical considerations for setting up QM/MM simulations of catalytic systems, are introduced; then, areas of heterogeneous catalysis are explored where QM/MM methods have been most fruitfully applied. The discussion includes simulations performed for adsorption processes in solvent at metallic interfaces, reaction mechanisms within zeolitic systems, nanoparticles, and defect chemistry within ionic solids. We conclude with a perspective on the current state of the field and areas where future opportunities for development and application exist.

Distal Mutations in the β-Clamp of DNA Polymerase III* Disrupt DNA Orientation and Affect Exonuclease Activity

DNA polymerases are responsible for the replication and repair of DNA found in all DNA-based organisms. DNA Polymerase III is the main replicative polymerase of E. coli and is composed of over 10 proteins. A subset of these proteins (Pol III*) includes the polymerase (α), exonuclease (ϵ), clamp (β), and accessory protein (θ). Mutations of residues in, or around the active site of the catalytic subunits (α and ϵ), can have a significant impact on catalysis. However, the effects of distal mutations in noncatalytic subunits on the activity of catalytic subunits are less well-characterized. Here, we investigate the effects of two Pol III* variants, β-L82E/L82'E and β-L82D/L82'D, on the proofreading reaction catalyzed by ϵ. MD simulations reveal major changes in the dynamics of Pol III*, which extend throughout the complex. These changes are mostly induced by a shift in the position of the DNA substrate inside the β-clamp, although no major structural changes are observed in the protein complex. Quantum mechanics/molecular mechanics (QM/MM) calculations indicate that the β-L82D/L82'D variant has reduced catalytic proficiency due to highly endoergic reaction energies resulting from structural changes in the active site and differences in the electric field at the active site arising from the protein and substrate. Conversely, the β-L82E/L82'E variant is predicted to maintain proofreading activity, exhibiting a similar reaction barrier for nucleotide excision compared with the WT system. However, significant differences in the reaction mechanism are obtained due to the changes induced by the mutations on the β-clamp.

Computational Investigation of a Series of Small Molecules as Potential Compounds for Lysyl Hydroxylase-2 (LH2) Inhibition

The catalytic function of lysyl hydroxylase-2 (LH2), a member of the Fe(II)/αKG-dependent oxygenase superfamily, is to catalyze the hydroxylation of lysine to hydroxylysine in collagen, resulting in stable hydroxylysine aldehyde-derived collagen cross-links (HLCCs). Reports show that high amounts of LH2 lead to the accumulation of HLCCs, causing fibrosis and specific types of cancer metastasis. Some members of the Fe(II)/αKG-dependent family have also been reported to have intramolecular O2 tunnels, which aid in transporting one of the required cosubstrates into the active site. While LH2 can be a promising target to combat these diseases, efficacious inhibitors are still lacking. We have used computational simulations to investigate a series of 44 small molecules as lead compounds for LH2 inhibition. Tunneling analyses indicate the existence of several intramolecular tunnels. The lengths of the calculated O2-transporting tunnels in holoenzymes are relatively longer than those in the apoenzyme, suggesting that the ligands may affect the enzyme's structure and possibly block (at least partially) the tunnels. The sequence alignment analysis between LH enzymes from different organisms shows that all of the amino acid residues with the highest occurrence rate in the oxygen tunnels are conserved. Our results suggest that the enolate form of diketone compounds establishes stronger interactions with the Fe(II) in the active site. Branching the enolate compounds with functional groups such as phenyl and pyridinyl enhances the interaction with various residues around the active site. Our results provide information about possible leads for further LH2 inhibition design and development.

Electronic Structure Effects Related to the Origin of the Remarkable Near-Infrared Absorption of Blastochloris viridis' Light Harvesting 1-Reaction Center Complex

Various photosynthetic organisms have evolved to absorb light in different regions of the visible light spectrum, thus adapting to the various lighting conditions available on Earth. While most of these autotrophic organisms absorb wavelengths around the 700-800 nm region, some are capable of red-shifted absorptions above this range, but none as remarkably as Blastochloris viridis whose main absorption is observed at 1015 nm, approximately 220 nm (0.34 eV) lower in energy than their main constituent pigments, BChl-b, whose main absorption is observed at 795 nm. The structure of its light harvesting 1-reaction center was recently elucidated by cryo-EM; however, the electronic structure details behind this red-shifted absorption remain unattended. We used hybrid quantum mechanics/molecular mechanics (QM/MM) calculations to optimize one of the active centers and performed classical molecular dynamics (MD) simulations to sample conformations beyond the optimized structure. We did excited state calculations with the time-dependent density functional theory method at the CAM-B3LYP/cc-pVDZ level of theory. We reproduced the near IR absorption by sequentially modifying the number of components involved in our systems using representative structures from the calculated MD ensemble. Natural transition orbital analysis reveals the participation of the BChl-b fragments to the main transition in the native structure and the structures obtained from the QM/MM and MD simulations. H-bonding pigment-protein interactions play a role on the conformation stabilization and orientation; however, the bacteriochlorin ring conformations and the exciton delocalization are the most relevant factors to explain the red-shifting phenomenon.

Atomistic-Level Description of the Covalent Inhibition of SARS-CoV-2 Papain-like Protease

Inhibition of the papain-like protease (PLpro) of SARS-CoV-2 has been demonstrated to be a successful target to prevent the spreading of the coronavirus in the infected body. In this regard, covalent inhibitors, such as the recently proposed VIR251 ligand, can irreversibly inactivate PLpro by forming a covalent bond with a specific residue of the catalytic site (Cys111), through a Michael addition reaction. An inhibition mechanism can therefore be proposed, including four steps: (i) ligand entry into the protease pocket; (ii) Cys111 deprotonation of the thiol group by a Brønsted-Lowry base; (iii) Cys111-S- addition to the ligand; and (iv) proton transfer from the protonated base to the covalently bound ligand. Evaluating the energetics and PLpro conformational changes at each of these steps could aid the design of more efficient and selective covalent inhibitors. For this aim, we have studied by means of MD simulations and QM/MM calculations the whole mechanism. Regarding the first step, we show that the inhibitor entry in the PLpro pocket is thermodynamically favorable only when considering the neutral Cys111, that is, prior to the Cys111 deprotonation. For the second step, MD simulations revealed that His272 would deprotonate Cys111 after overcoming an energy barrier of ca. 32 kcal/mol (at the QM/MM level), but implying a decrease of the inhibitor stability inside the protease pocket. This information points to a reversible Cys111 deprotonation, whose equilibrium is largely shifted toward the neutral Cys111 form. Although thermodynamically disfavored, if Cys111 is deprotonated in close proximity to the vinylic carbon of the ligand, then covalent binding takes place in an irreversible way (third step) to form the enolate intermediate. Finally, due to Cys111-S- negative charge redistribution over the bound ligand, proton transfer from the initially protonated His272 is favored, finally leading to an irreversibly modified Cys111 and a restored His272. These results elucidate the selectivity of Cys111 to enable formation of a covalent bond, even if a weak proton acceptor is available, as His272.

Computational investigations of selected enzymes from two iron and α-ketoglutarate-dependent families

DNA alkylation is used as the key epigenetic mark in eukaryotes, however, most alkylation in DNA can result in deleterious effects. Therefore, this process needs to be tightly regulated. The enzymes of the AlkB and Ten-Eleven Translocation (TET) families are members of the Fe and alpha-ketoglutarate-dependent superfamily of enzymes that are tasked with dealkylating DNA and RNA in cells. Members of these families span all species and are an integral part of transcriptional regulation. While both families catalyze oxidative dealkylation of various bases, each has specific preference for alkylated base type as well as distinct catalytic mechanisms. This perspective aims to provide an overview of computational work carried out to investigate several members of these enzyme families including AlkB, ALKB Homolog 2, ALKB Homolog 3 and Ten-Eleven Translocate 2. Insights into structural details, mutagenesis studies, reaction path analysis, electronic structure features in the active site, and substrate preferences are presented and discussed.

Structural and electronic analysis of the octarepeat region of prion protein with four Cu2+ by polarizable MD and QM/MM simulations

Prions have been linked to neurodegenerative diseases that affect various species of mammals including humans. The prion protein, located mainly in neurons, is believed to play the role of metal ion transporter. High levels of copper ions have been related to structural changes. A 32-residue region of the N-terminal domain, known as octarepeat, can bind up to four copper ions. Different coordination modes have been observed and are strongly dependent on Cu2+ concentration. Many theoretical studies carried out so far have focused on studying the coordination modes of a single copper ion. In this work we investigate the octarepeat region coordinated with four copper ions. Molecular dynamics (MD) and hybrid quantum mechanics/molecular mechanics (QM/MM) simulations using the polarizable AMOEBA force field have been carried out. The polarizable MD simulations starting from a fully extended conformation indicate that the tetra-Cu2+/octarepeat complex forms a globular structure. The globular form is stabilized by interactions between Cu2+ and tryptophan residues resulting in some coordination sites observed to be in close proximity, in agreement with experimental results. Subsequent QM/MM simulations on several snapshots suggests the system is in a high-spin quintet state, with all Cu2+ bearing one single electron, and all unpaired electrons are ferromagnetically coupled. NMR simulations on selected structures provides insights on the chemical shifts of the first shell ligands around the metals with respect to inter-metal distances.

Combining Evolutionary Conservation and Quantum Topological Analyses To Determine Quantum Mechanics Subsystems for Biomolecular Quantum Mechanics/Molecular Mechanics Simulations

Selection of residues and other molecular fragments for inclusion in the quantum mechanics (QM) region for QM/molecular mechanics (MM) simulations is an important step for these calculations. Here, we present an approach that combines protein sequence/structure evolution and electron localization function (ELF) analyses. The combination of these two analyses allows the determination of whether a residue needs to be included in the QM subsystem or can be represented by the MM environment. We have applied this approach on two systems previously investigated by QM/MM simulations, 4-oxalocrotonate tautomerase (4OT) and ten-eleven translocation-2 (TET2), that provide examples where fragments may or may not need to be included in the QM subsystem. Subsequently, we present the use of this approach to determine the appropriate QM subsystem to calculate the minimum energy path (MEP) for the reaction catalyzed by human DNA polymerase λ (Polλ) with a third cation in the active site. Our results suggest that the combination of protein evolutionary and ELF analyses provides insights into residue/molecular fragment selection for QM/MM simulations.

Development and application of quantum mechanics/molecular mechanics methods with advanced polarizable potentials

Quantum mechanics/molecular mechanics (QM/MM) simulations are a popular approach to study various features of large systems. A common application of QM/MM calculations is in the investigation of reaction mechanisms in condensed-phase and biological systems. The combination of QM and MM methods to represent a system gives rise to several challenges that need to be addressed. The increase in computational speed has allowed the expanded use of more complicated and accurate methods for both QM and MM simulations. Here, we review some approaches that address several common challenges encountered in QM/MM simulations with advanced polarizable potentials, from methods to account for boundary across covalent bonds and long-range effects, to polarization and advanced embedding potentials.

A Many-Body, Fully Polarizable Approach to QM/MM Simulations

We present a new development in quantum mechanics/molecular mechanics (QM/MM) methods by replacing conventional MM models with data-driven many-body (MB) representations rigorously derived from high-level QM calculations. The new QM/MM approach builds on top of mutually polarizable QM/MM schemes developed for polarizable force fields with inducible dipoles and uses permutationally invariant polynomials to effectively account for quantum-mechanical contributions (e.g., exchange-repulsion and charge transfer and penetration) that are difficult to describe by classical expressions adopted by conventional MM models. Using the many-body MB-pol and MB-DFT potential energy functions for water, which include explicit two-body and three-body terms fitted to reproduce the corresponding CCSD(T) and PBE0 two-body and three-body energies for water, we demonstrate a smooth energetic transition as molecules are transferred between QM and MM regions, without the need of a transition layer. By effectively elevating the accuracy of both the MM region and the QM/MM interface to that of the QM region, the new QM/MB-MM approach achieves an accuracy comparable to that obtained with a fully QM treatment of the entire system.

Ewald-based methods for Gaussian integral evaluation: application to a new parameterization of GEM

The development of accurate potentials for computational simulations has been an active area of research. Our group has been involved in the development of the Gaussian electrostatic model (GEM), a force field based on molecular densities. The philosophy of GEM is based on the pioneering work of N. Gresh and co-workers of the reproduction of individual inter-molecular interaction components obtained from quantum mechanical (QM) energy decomposition analysis (EDA). The molecular densities used in GEM are represented by fitting accurate QM molecular densities using auxiliary basis sets (comprised of Hermite Gaussians). The use of these molecular densities results in the need to evaluate a large number of Gaussian integrals. We have previously shown that the particle-mesh Ewald (PME), and fast Fourier Poisson (FFP) methods can be used for efficiently evaluating these types of integrals. Here, we present the latest parameterization of GEM* and its application for an extensive study of PME and FFP for molecular dynamics (MD) simulations using a hybrid version of our potential, GEM*. The temperature dependence of various bulk properties is presented and discussed, as well as the effect of various parameters affecting the performance/accuracy of both methods.

Ground State Destabilization in Uracil DNA Glycosylase: Let's Not Forget "Tautomeric Strain" in Substrates

Enzymes like uracil DNA glycosylase (UDG) can achieve ground state destabilization, by polarizing substrates to mimic rare tautomers. On the basis of computed nucleus independent chemical shifts, NICS(1)_zz, and harmonic oscillator model of electron delocalization (HOMED) analyses, of quantum mechanics (QM) and quantum mechanics/molecular mechanics (QM/MM) models of the UDG active site, uracil is strongly polarized when bound to UDG and resembles a tautomer >12 kcal/mol higher in energy. Natural resonance theory (NRT) analyses identified a dominant O2 imidate resonance form for residue bound 1-methyl-uracil. This "tautomeric strain" raises the energy of uracil, making uracilate a better than expected leaving group. Computed gas-phase S_N2 reactions of free and hydrogen bonded 1-methyl-uracil demonstrate the relationship between the degree of polarization in uracil and the leaving group ability of uracilate.

Towards large scale hybrid QM/MM dynamics of complex systems with advanced point dipole polarizable embeddings

In this work, we present a general route to hybrid Quantum Mechanics/Molecular Mechanics (QM/MM) Molecular Dynamics for complex systems using a polarizable embedding. We extend the capabilities of our hybrid framework, combining the Gaussian and Tinker/Tinker-HP packages in the context of the AMOEBA polarizable force field to treat large (bio)systems where the QM and the MM subsystems are covalently bound, adopting pseudopotentials at the boundaries between the two regions. We discuss in detail the implementation and demonstrate the global energy conservation of our QM/MM Born-Oppenheimer molecular dynamics approach using Density Functional Theory. Finally, the approach is assessed on the electronic absorption properties of a 16 500 atom complex encompassing an organic dye embedded in a DNA matrix in solution, extending the hybrid method to a time-dependent Density Functional Theory approach. The results obtained comparing different partitions between the quantum and the classical subsystems also suggest that large QM portions are not necessary if accurate polarizable force fields are used in a variational formulation of the embedding, properly including the QM/MM mutual polarization.