Exploring the Dependence of QM/MM Calculations of Enzyme Catalysis on the Size of the QM Region

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.

Condensed Matter Systems Exposed to Radiation: Multiscale Theory, Simulations, and Experiment

This roadmap reviews the new, highly interdisciplinary research field studying the behavior of condensed matter systems exposed to radiation. The Review highlights several recent advances in the field and provides a roadmap for the development of the field over the next decade. Condensed matter systems exposed to radiation can be inorganic, organic, or biological, finite or infinite, composed of different molecular species or materials, exist in different phases, and operate under different thermodynamic conditions. Many of the key phenomena related to the behavior of irradiated systems are very similar and can be understood based on the same fundamental theoretical principles and computational approaches. The multiscale nature of such phenomena requires the quantitative description of the radiation-induced effects occurring at different spatial and temporal scales, ranging from the atomic to the macroscopic, and the interlinks between such descriptions. The multiscale nature of the effects and the similarity of their manifestation in systems of different origins necessarily bring together different disciplines, such as physics, chemistry, biology, materials science, nanoscience, and biomedical research, demonstrating the numerous interlinks and commonalities between them. This research field is highly relevant to many novel and emerging technologies and medical applications.

CHARMM-GUI QM/MM Interfacer for a Quantum Mechanical and Molecular Mechanical (QM/MM) Simulation Setup: 1. Semiempirical Methods

Quantum mechanical (QM) treatments, when combined with molecular mechanical (MM) force fields, can effectively handle enzyme-catalyzed reactions without significantly increasing the computational cost. In this context, we present CHARMM-GUI QM/MM Interfacer, a web-based cyberinfrastructure designed to streamline the preparation of various QM/MM simulation inputs with ligand modification. The development of QM/MM Interfacer has been achieved through integration with existing CHARMM-GUI modules, such as PDB Reader and Manipulator, Solution Builder, and Membrane Builder. In addition, new functionalities have been developed to facilitate the one-stop preparation of QM/MM systems and enable interactive and intuitive ligand modifications and QM atom selections. QM/MM Interfacer offers support for a range of semiempirical QM methods, including AM1(+/d), PM3(+/PDDG), MNDO(+/d, +/PDDG), PM6, RM1, and SCC-DFTB, tailored for both AMBER and CHARMM. A nontrivial setup related to ligand modification, link-atom insertion, and charge distribution is automatized through intuitive user interfaces. To illustrate the robustness of QM/MM Interfacer, we conducted QM/MM simulations of three enzyme-substrate systems: dihydrofolate reductase, insulin receptor kinase, and oligosaccharyltransferase. In addition, we have created three tutorial videos about building these systems, which can be found at https://www.charmm-gui.org/demo/qmi. QM/MM Interfacer is expected to be a valuable and accessible web-based tool that simplifies and accelerates the setup process for hybrid QM/MM simulations.

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.

Catalytic Effects of Active Site Conformational Change in the Allosteric Activation of Imidazole Glycerol Phosphate Synthase

Imidazole glycerol phosphate synthase (IGPS) is a class-I glutamine amidotransferase (GAT) that hydrolyzes glutamine. Ammonia is produced and transferred to a second active site, where it reacts with N1-(5'-phosphoribosyl)-formimino-5-aminoimidazole-4-carboxamide ribonucleotide (PrFAR) to form precursors to purine and histidine biosynthesis. Binding of PrFAR over 25 Å away from the active site increases glutaminase efficiency by ∼4500-fold, primarily altering the glutamine turnover number. IGPS has been the focus of many studies on allosteric communication; however, atomic details for how the glutamine hydrolysis rate increases in the presence of PrFAR are lacking. We present a density functional theory study on 237-atom active site cluster models of IGPS based on crystallized structures representing the inactive and allosterically active conformations and investigate the multistep reaction leading to thioester formation and ammonia production. The proposed mechanism is supported by similar, well-studied enzyme mechanisms, and the corresponding energy profile is consistent with steady-state kinetic studies of PrFAR + IGPS. Additional active site models are constructed to examine the relationship between active site structural change and transition-state stabilization via energy decomposition schemes. The results reveal that the inactive IGPS conformation does not provide an adequately formed oxyanion hole structure and that repositioning of the oxyanion strand relative to the substrate is vital for a catalysis-competent oxyanion hole, with or without the hVal51 dihedral flip. These findings are valuable for future endeavors in modeling the IGPS allosteric mechanism by providing insight into the atomistic changes required for rate enhancement that can inform suitable reaction coordinates for subsequent investigations.

Peptidomimetics for CVD screened via TRADD-TRAF2 complex interface assessments

The main aim of this study is to screen and develop Peptidomimetics to treat atherosclerosis (AS) which is a Cardio Vascular Disease (CVD). Peptidomimetics were obtained from the protein-protein interaction interface of TRADD (Tumor necrosis factor receptor type 1-associated DEATH domain protein) and TRAF2 (TNF receptor-associated factor 2) complex. TRADD-TRAF2 interaction is critical in AS pathogenesis since it assists a series of signal transducers that activate NF-κB. Conceptually, the triggered NF-κB makes an extensive amount of nitric oxide (NO) synthesized by inducible nitric oxide synthase (iNOS), which boons the progress of AS. The examined TRADD-TRAF2 complex (PDB ID: 1F3V) and its interaction details revealed that the sequence range W11-G165 of TRADD highly interacts with TRAF2. The sequence range W11-G165 was selected for the design and preparation of the inhibitory peptide in silico. The selected sequence was mutated with the alanine scanning method to have a range of inhibitory peptides. With the help of different in silico tools, the top three, namely MIP11-25 L, MIP131-143 h, and MIP149-164 m peptides showed the best interaction with the critical residues of TRAF2. Thus, these three peptides were used for generating peptidomimetics using pepMMsMIMIC, a peptidomimetics virtual screening tool. Around 600 peptidomimetics were identified & and retrieved for further screening by employing molecular docking tools and MD (Molecular Dynamics) simulations. Density Functional Theory (DFT) and ADMET predictions were applied to validate the screened peptidomimetics druggability. In the results, peptidomimic compounds MMs03918858 and MMs03927281 with binding energy values of -9.6 kcal/mol and - 9.1 kcal/mol respectively were screened as the best and are proposed for further pre-clinical studies.

Unraveling the Pivotal Roles of Various Metal Ion Centers in the Catalysis of Quercetin 2,4-Dioxygenases

Quercetin 2,4-dioxygenase (QueD) with various transition metal ion co-factors shows great differences, but the internal reasons have not been illustrated in detail. In order to explore the effects of metal ion centers on the catalytic reactivity of QueD, we calculated and compared the minimum energy crossing point (MECP) of dioxygen from the relatively stable triplet state to the active singlet state under different conditions by using the DFT method. It was found that the metal ions play a more important role in the activation of dioxygen compared with the substrate and the protein environment. Simultaneously, the catalytic reactions of the bacterial QueDs containing six different transition metal ions were studied by the QM/MM approach, and we finally obtained the reactivity sequence of metal ions, Ni2+ > Co2+ > Zn2+ > Mn2+ > Fe2+ > Cu2+, which is basically consistent with the previous experimental results. Our calculation results indicate that metal ions act as Lewis acids in the reaction to stabilize the substrate anion and the subsequent superoxo and peroxo species in the reaction, and promote the proton coupled electron transfer (PCET) process. Furthermore, the coordination tendencies of transition metal ion centers also have important effects on the catalytic cycle. These findings have general implications on metalloenzymes, which can expand our understanding on how various metal ions play their key role in modulating catalytic reactivity.

High-throughput computational investigation of protein electrostatics and cavity for SAM-dependent methyltransferases

S-adenosyl methionine (SAM)-dependent methyl transferases (MTases) are a ubiquitous class of enzymes catalyzing dozens of essential life processes. Despite targeting a large space of substrates with diverse intrinsic reactivity, SAM MTases have similar catalytic efficiency. While understanding of MTase mechanism has grown tremendously through the integration of structural characterization, kinetic assays, and multiscale simulations, it remains elusive how these enzymes have evolved to fit the diverse chemical needs of their respective substrates. In this work, we performed a high-throughput molecular modeling analysis of 91 SAM MTases to better understand how their properties (i.e., electric field [EF] strength and active site volumes) help achieve similar catalytic efficiency toward substrates of different reactivity. We found that EF strengths have largely adjusted to make the target atom a better methyl acceptor. For MTases that target RNA/DNA and histone proteins, our results suggest that EF strength accommodates formal hybridization state and variation in cavity volume trends with diversity of substrate classes. Metal ions in SAM MTases contribute negatively to EF strength for methyl donation and enzyme scaffolds tend to offset these contributions.

Efficient automatic construction of atom-economical QM regions with point-charge variation analysis

The setup of QM/MM calculations is not trivial since many decisions have to be made by the simulation scientist to achieve reasonable and consistent results. The main challenge to be tackled is the construction of the QM region to make sure to take into account all important parts of the adjacent environment and exclude less important ones. In our previous work [F. Brandt and Ch. R. Jacob, Systematic QM Region Construction in QM/MM Calculations Based on Uncertainty Quantification, J. Chem. Theory Comput., 2022, 18, 2584-2596.], we introduced the point charge variation analysis (PCVA) as a simple and reliable tool to systematically construct QM regions based on the sensitivity of the reaction energy with respect to variations of the MM point charges. Here, we assess several simplified variants of this PCVA approach for the example of catechol O-methyltransferase and apply PCVA for another system, the triosephosphate isomerase. Furthermore, we extend its scope by applying it to a DNA system. Our results indicate that PCVA offers an efficient and versatile approach of the automatic construction of atom-economical QM regions, but also identify possible pitfalls and limitations.

Hybrid Quantum Mechanical/Molecular Mechanical Methods For Studying Energy Transduction in Biomolecular Machines

Hybrid quantum mechanical/molecular mechanical (QM/MM) methods have become indispensable tools for the study of biomolecules. In this article, we briefly review the basic methodological details of QM/MM approaches and discuss their applications to various energy transduction problems in biomolecular machines, such as long-range proton transports, fast electron transfers, and mechanochemical coupling. We highlight the particular importance for these applications of balancing computational efficiency and accuracy. Using several recent examples, we illustrate the value and limitations of QM/MM methodologies for both ground and excited states, as well as strategies for calibrating them in specific applications. We conclude with brief comments on several areas that can benefit from further efforts to make QM/MM analyses more quantitative and applicable to increasingly complex biological problems.

Fragment-Based Calculations of Enzymatic Thermochemistry Require Dielectric Boundary Conditions

Electronic structure calculations on enzymes require hundreds of atoms to obtain converged results, but fragment-based approximations offer a cost-effective solution. We present calculations on enzyme models containing 500-600 atoms using the many-body expansion, comparing to benchmarks in which the entire enzyme-substrate complex is described at the same level of density functional theory. When the amino acid fragments contain ionic side chains, the many-body expansion oscillates under vacuum boundary conditions but rapid convergence is restored using low-dielectric boundary conditions. This implies that full-system calculations in the gas phase are inappropriate benchmarks for assessing errors in fragment-based approximations. A three-body protocol retains sub-kilocalorie per mole fidelity with respect to a supersystem calculation, as does a two-body calculation combined with a full-system correction at a low-cost level of theory. These protocols pave the way for application of high-level quantum chemistry to large systems via rigorous, ab initio treatment of many-body polarization.

A Multiple Proton Transfer Mechanism for the Charging Step of the Aminoacylation Reaction at the Active Site of Aspartyl tRNA Synthetase

Aspartyl-tRNA synthetase catalyzes the attachment of aspartic acid to its cognate tRNA by the aminoacylation reaction during the initiation of the protein biosynthesis process. In the second step of the aminoacylation reaction, known as the charging step, the aspartate moiety is transferred from aspartyl-adenylate to the 3'-OH of A76 of tRNA through a proton transfer process. We have investigated different pathways for the charging step through three separate QM/MM simulations combined with the enhanced sampling method of well-sliced metadynamics and found out the most feasible pathway for the reaction at the active site of the enzyme. In the charging reaction, both the phosphate group and the ammonium group after deprotonation can potentially act as a base for proton transfer in the substrate-assisted mechanism. We have considered three possible mechanisms involving different pathways of proton transfer, and only one of them is determined to be enzymatically feasible. The free energy landscape along reaction coordinates where the phosphate group acts as the general base showed that, in the absence of water, the barrier height is 52.6 kcal/mol. The free energy barrier is reduced to 39.7 kcal/mol when the active site water molecules are also treated quantum mechanically, thus allowing a water mediated proton transfer. The charging reaction involving the ammonium group of the aspartyl adenylate is found to follow a path where first a proton from the ammonium group moves to a water in the vicinity forming a hydronium ion (H₃O⁺) and NH₂ group. The hydronium ion subsequently passes the proton to the Asp233 residue, thus minimizing the chance of back proton transfer from hydronium to the NH₂ group. The neutral NH₂ group subsequently takes the proton from the O3' of A76 with a free energy barrier of 10.7 kcal/mol. In the next step, the deprotonated O3' makes a nucleophilic attack to the carbonyl carbon forming a tetrahedral transition state with a free energy barrier of 24.8 kcal/mol. Thus, the present work shows that the charging step proceeds through a multiple proton transfer mechanism where the amino group formed after deprotonation acts as the base to capture a proton from O3' of A76 rather than the phosphate group. The current study also shows the important role played by Asp233 in the proton transfer process.

Accelerated Quantum Mechanics/Molecular Mechanics Simulations via Neural Networks Incorporated with Mechanical Embedding Scheme

A powerful tool to study the mechanism of reactions in solutions or enzymes is to perform the ab initio quantum mechanical/molecular mechanical (QM/MM) molecular dynamics (MD) simulations. However, the computational cost is too high due to the explicit electronic structure calculations at every time step of the simulation. A neural network (NN) method can accelerate the QM/MM-MD simulations, but it has long been a problem to accurately describe the QM/MM electrostatic coupling by NN in the electrostatic embedding (EE) scheme. In this work, we developed a new method to accelerate QM/MM calculations in the mechanic embedding (ME) scheme. The potentials and partial point charges of QM atoms are first learned in vacuo by the embedded atom neural networks (EANN) approach. MD simulations are then performed on this EANN/MM potential energy surface (PES) to obtain free energy (FE) profiles for reactions, in which the QM/MM electrostatic coupling is treated in the mechanic embedding (ME) scheme. Finally, a weighted thermodynamic perturbation (wTP) corrects the FE profiles in the ME scheme to the EE scheme. For two reactions in water and one in methanol, our simulations reproduced the B3LYP/MM free energy profiles within 0.5 kcal/mol with a speed-up of 30-60-fold. The results show that the strategy of combining EANN potential in the ME scheme with the wTP correction is efficient and reliable for chemical reaction simulations in liquid. Another advantage of our method is that the QM PES is independent of the MM subsystem, so it can be applied to various MM environments as demonstrated by an S_N2 reaction studied in water and methanol individually, which used the same EANN PES. The free energy profiles are in excellent accordance with the results obtained from B3LYP/MM-MD simulations. In future, this method will be applied to the reactions of enzymes and their variants.

Insights into the Degradation of Polymer-Drug Conjugates by an Overexpressed Enzyme in Cancer Cells

Intensive efforts have been made to provide better treatments to cancer patients. Currently, nanoparticle-based drug delivery systems have gained propulsion, as they can overcome the drawbacks of free drugs. However, drug stability inside the nanocapsule must be ensured to prevent burst release. To overcome this, drugs conjugated to amphiphilic copolymers, assembled into nanoparticles, can provide a sustained release if endogenously degraded. Thus, we have designed and assessed the drug release viability of polymer-drug conjugates by the human Carboxylesterase 2, for a targeted drug activation. We performed molecular dynamics simulations applying a quantum mechanics/molecular mechanics potential to study the degradation profiles of 30 designed conjugates, where six were predicted to be hydrolyzed by this enzyme. We further analyzed the enzyme-substrate environment to delve into what structural features may lead to successful hydrolysis. These findings contribute to the development of new medicines ensuring effective cancer treatments with fewer side effects.

Challenges in Protein QM/MM Simulations with Intra-Backbone Link Atoms

Hybrid quantum mechanical/molecular mechanical (QM/MM) simulations fuel discoveries in many fields of science including computational biochemistry and enzymology. Development of more convenient tools leads to an increase in the number of works in which mechanical insights into enzymes' mode of operation are obtained. Most commonly, these tools feature hydrogen-capping (link atom) approach to provide coupling between QM and MM subsystems across a covalent bond. Extensive studies were conducted to provide a solid foundation for the correctness of such an approach when a bond to a nonpolar MM atom is considered. However, not every task may be accomplished this way. Certain scenarios of using QM/MM in computational enzymology encourage or even necessitate the incorporation of backbone atoms into the QM region. Two out of three backbone atoms are polar, and in QM/MM with electrostatic embedding, a neighboring link atom will be hyperpolarized. Several schemes to mitigate this effect were previously proposed alongside a rigorous assessment of quantitative effects on model systems. However, it was not clear whether they may translate into qualitatively different results and how link atom hyperpolarization may manifest itself in a real-life enzymological scenario. Here, we show that the consequences of such an artifact may be severe and may completely overturn the conclusions drawn from the simulations. Our case advocates for the use of charge redistribution schemes whenever intra-backbone QM/MM boundaries are considered. Moreover, we addressed how different boundary types and charge redistribution schemes influence backbone dynamics. We showed that the results are heavily dependent on which boundary MM terms are retained, with charge alteration being of secondary importance. In the worst case, only three intra-backbone boundaries may be used with relative confidence in the adequacy of resulting simulations, irrespective of the hyperpolarization mitigation scheme. Thus, advances in the field are certainly needed to fuel new discoveries. As of now, we believe that issues raised in this work might encourage authors in the field to report what boundaries, boundary MM terms, and charge redistribution schemes they are using, so their results may be correctly interpreted.

Effect of the QM Size, Basis Set, and Polarization on QM/MM Interaction Energy Decomposition Analysis

Herein, an Energy Decomposition Analysis (EDA) scheme extended to the framework of QM/MM calculations in the context of electrostatic embeddings (QM/MM-EDA) including atomic charges and dipoles is applied to assess the effect of the QM region size on the convergence of the different interaction energy components, namely, electrostatic, Pauli, and polarization, for cationic, anionic, and neutral systems interacting with a strong polar environment (water). Significant improvements are found when the bulk solvent environment is described by a MM potential in the EDA scheme as compared to pure QM calculations that neglect bulk solvation. The predominant electrostatic interaction requires sizable QM regions. The results reported here show that it is necessary to include a surprisingly large number of water molecules in the QM region to obtain converged values for this energy term, contrary to most cluster models often employed in the literature. Both the improvement of the QM wave function by means of a larger basis set and the introduction of polarization into the MM region through a polarizable force field do not translate to a faster convergence with the QM region size, but they lead to better results for the different interaction energy components. The results obtained in this work provide insight into the effect of each energy component on the convergence of the solute-solvent interaction energy with the QM region size. This information can be used to improve the MM FFs and embedding schemes employed in QM/MM calculations of solvated systems.

Convergence in determining enzyme functional descriptors across Kemp eliminase variants

Molecular simulations have been extensively employed to accelerate biocatalytic discoveries. Enzyme functional descriptors derived from molecular simulations have been leveraged to guide the search for beneficial enzyme mutants. However, the ideal active-site region size for computing the descriptors over multiple enzyme variants remains untested. Here, we conducted convergence tests for dynamics-derived and electrostatic descriptors on 18 Kemp eliminase variants across six active-site regions with various boundary distances to the substrate. The tested descriptors include the root-mean-square deviation of the active-site region, the solvent accessible surface area ratio between the substrate and active site, and the projection of the electric field (EF) on the breaking C-H bond. All descriptors were evaluated using molecular mechanics methods. To understand the effects of electronic structure, the EF was also evaluated using quantum mechanics/molecular mechanics methods. The descriptor values were computed for 18 Kemp eliminase variants. Spearman correlation matrices were used to determine the region size condition under which further expansion of the region boundary does not substantially change the ranking of descriptor values. We observed that protein dynamics-derived descriptors, including RMSD_{active_site} and SASA_ratio, converge at a distance cutoff of 5 Å from the substrate. The electrostatic descriptor, EF_C-H, converges at 6 Å using molecular mechanics methods with truncated enzyme models and 4 Å using quantum mechanics/molecular mechanics methods with whole enzyme model. This study serves as a future reference to determine descriptors for predictive modeling of enzyme engineering.

Assessment of Reversibility for Covalent Cysteine Protease Inhibitors Using Quantum Mechanics/Molecular Mechanics Free Energy Surfaces

We have used molecular dynamics (MD) simulations with hybrid quantum mechanics/molecular mechanics (QM/MM) potentials to investigate the reaction mechanism for covalent inhibition of cathepsin K and assess the reversibility of inhibition. The computed free energy profiles suggest that a nucleophilic attack by the catalytic cysteine on the inhibitor warhead and proton transfer from the catalytic histidine occur in a concerted manner. The results indicate that the reaction is more strongly exergonic for the alkyne-based inhibitors, which bind irreversibly to cathepsin K, than for the nitrile-based inhibitor odanacatib, which binds reversibly. Gas-phase energies were also calculated for the addition of methanethiol to structural prototypes for a number of warheads of interest in cysteine protease inhibitor design in order to assess electrophilicity. The approaches presented in this study are particularly applicable to assessment of novel warheads, and computed transition state geometries can be incorporated into molecular models for covalent docking.

Improving the Accuracy of Quantum Mechanics/Molecular Mechanics (QM/MM) Models with Polarized Fragment Charges

This paper introduces an economical approach for improving the accuracy and convergence of quantum mechanics/molecular mechanics (QM/MM) models. The approach is tested on a series of neutral and charged amino acids embedded in a 160-water cluster, where their intramolecular proton transfer energies (neutral amino acid → zwitterionic amino acid) were previously obtained at the ωB97X-D/6-31G(d) level of theory. When the charges on the MM atoms were replaced with those obtained at the same QM level of theory used to treat the QM atoms, this significantly improved the accuracy and convergence of the QM/MM models. In particular, the QM/MM model converged to within 1.4 kcal mol^-1 of directly calculated DFT energies for smaller (by as many as 20 waters) QM regions. The use of atomic charges obtained from the natural population analysis yielded the most significant improvement, while other charge schemes such as Mulliken, electrostatic potential, or CM5 led to poorer outcomes. It is further demonstrated that the QM atomic charges can be accurately estimated in a highly efficient manner using an iterative fragmentation approach based on the moving-domain QM/MM method. Similar observations were made when the approach was used to predict the barrier of an S_N2 reaction. Thus, the use of QM-quality atomic charges on MM atoms represents a simple and easy-to-implement strategy for improving the accuracy of QM/MM models.

Free Energy Landscape of the Adenylation Reaction of the Aminoacylation Process at the Active Site of Aspartyl tRNA Synthetase

The process of protein biosynthesis is initiated by the aminoacylation process where a transfer ribonucleic acid (tRNA) is charged by the attachment of its cognate amino acid at the active site of the corresponding aminoacyl tRNA synthetase enzyme. The first step of the aminoacylation process, known as the adenylation reaction, involves activation of the cognate amino acid where it reacts with a molecule of adenosine triphosphate (ATP) at the active site of the enzyme to form the aminoacyl adenylate and inorganic pyrophosphate. In the current work, we have investigated the adenylation reaction between aspartic acid and ATP at the active site of the fully solvated aspartyl tRNA synthetase (AspRS) from Escherichia coli in aqueous medium at room temperature through hybrid quantum mechanical/molecular mechanical (QM/MM) simulations combined with enhanced sampling methods of well-tempered and well-sliced metadynamics. The objective of the present work is to study the associated free energy landscape and reaction barrier and also to explore the effects of active site mutation on the free energy surface of the reaction. The current calculations include finite temperature effects on free energy profiles. In particular, apart from contributions of interaction energies, the current calculations also include contributions of conformational, vibrational, and translational entropy of active site residues, substrates, and also the rest of the solvated protein and surrounding water into the free energy calculations. The present QM/MM metadynamics simulations predict a free energy barrier of 23.35 and 23.5 kcal/mol for two different metadynamics methods used to perform the reaction at the active site of the wild type enzyme. The free energy barrier increases to 30.6 kcal/mol when Arg217, which is an important conserved residue of the wild type enzyme at its active site, is mutated by alanine. These free energy results including the effect of mutation compare reasonably well with those of kinetic experiments that are available in the literature. The current work also provides molecular details of structural changes of the reactants and surroundings as the system dynamically evolves along the reaction pathway from reactant to the product state through QM/MM metadynamics simulations.

A Computational Study of the Promiscuity of the SAM-Dependent Methyltransferase AtHTMT1

A two-pronged computational approach was taken to study the promiscuity of the SAM⁺-dependent methyl transferase AtHTMT1 from thale cress with several nucleophiles (Cl^-, Br^-, I^-, NCO^-, NCS^-). First, enzyme-free methyl transfer reactions were studied with M05/6-311+G(2d,p) DFT calculations and electrostatic continuum models (PCM/SMD) for various chemical environments. Second, QM/MM MD simulations with semiempirical Hamiltonians (PM7, PM6-D3, AM1, PM6-D3H4) and the AMBER 14SB force field were used to study the enzyme catalyzed reaction in silico. The combination of the DFT and MD results shows that reactant desolvation generally accelerates the reaction, but it cannot explain the selectivity of the enzyme. The critical position of H₂O molecules at the reactive site favors the reaction of NCS^- over Cl^- and Br^- in agreement with experiments, but not observed in the quantum calculations for the cytosol. The addition of selected H₂O molecules to the N terminus of NCS^- greatly increases its reactivity, while H₂O molecules attached to Cl^- slow the reaction. The partial solvation of the nucleophiles in the reactive pouch holds the key to understanding the reactivity of AtHTMT1.

Systematic QM Region Construction in QM/MM Calculations Based on Uncertainty Quantification

While QM/MM studies of enzymatic reactions are widely used in computational chemistry, the results of such studies are subject to numerous sources of uncertainty, and the effect of different choices by the simulation scientist that are required when setting up QM/MM calculations is often unclear. In particular, the selection of the QM region is crucial for obtaining accurate and reliable results. Simply including amino acids by their distance to the active site is mostly not sufficient as necessary residues are missing or unimportant residues are included without evidence. Here, we take a first step toward quantifying uncertainties in QM/MM calculations by assessing the sensitivity of QM/MM reaction energies with respect to variations of the MM point charges. We show that such a point charge variation analysis (PCVA) can be employed to judge the accuracy of QM/MM reaction energies obtained with a selected QM region and devise a protocol to systematically construct QM regions that minimize this uncertainty. We apply such a PCVA to the example of catechol O-methyltransferase and demonstrate that it provides a simple and reliable approach for the construction of the QM region. Our PCVA-based scheme is computationally efficient and requires only calculations for a system with a minimal QM region. Our work highlights the promise of applying methods of uncertainty quantification in computational chemistry.

Factors That Determine the Variation of Equilibrium and Kinetic Properties of QM/MM Enzyme Simulations: QM Region, Conformation, and Boundary Condition

To analyze the impact of various technical details on the results of quantum mechanical (QM)/molecular mechanical (MM) enzyme simulations, including the QM region size, catechol-O-methyltransferase (COMT) is studied as a model system using an approximate QM/MM method (DFTB3/CHARMM). The results show that key equilibrium and kinetic properties for methyl transfer in COMT exhibit limited variations with respect to the size of the QM region, which ranges from ∼100 to ∼500 atoms in this study. With extensive sampling, local and global structural characteristics of the enzyme are largely conserved across the studied QM regions, while the nature of the transition state (e.g., secondary kinetic isotope effect) and reaction exergonicity are largely maintained. Deviations in the free energy profile with different QM region sizes are similar in magnitude to those observed with changes in other simulation protocols, such as different initial enzyme conformations and boundary conditions. Electronic structural properties, such as the covariance matrix of residual charge fluctuations, appear to exhibit rather long-range correlations, especially when the peptide backbone is included in the QM region; this observation holds when a range-separated DFT approach is used as the QM region, suggesting that delocalization error is unlikely the origin. Overall, the analyses suggest that multiple simulation details determine the results of QM/MM enzyme simulations with comparable contributions.

Hybrid QM/MM Free-Energy Evaluation of Drug-Resistant Mutational Effect on the Binding of an Inhibitor Indinavir to HIV-1 Protease

A human immunodeficiency virus-1 (HIV-1) protease is a homodimeric aspartic protease essential for the replication of HIV. The HIV-1 protease is a target protein in drug discovery for antiretroviral therapy, and various inhibitor molecules of transition state analogues have been developed. However, serious drug-resistant mutants have emerged. For understanding the molecular mechanism of the drug resistance, an accurate examination of the impacts of the mutations on ligand binding and enzymatic activity is necessary. Here, we present a molecular simulation study on the ligand binding of indinavir, a potent transition state analogue inhibitor, to the wild-type protein and a V82T/I84V drug-resistant mutant of the HIV-1 protease. We employed a hybrid ab initio quantum mechanical/molecular mechanical (QM/MM) free-energy optimization technique which combines a highly accurate QM description of the ligand molecule and its interaction with statistically ample conformational sampling of the MM protein environment by long-time molecular dynamics simulations. Through the free-energy calculations of protonation states of catalytic groups at the binding pocket and of the ligand-binding affinity changes upon the mutations, we successfully reproduced the experimentally observed significant reduction of the binding affinity upon the drug-resistant mutations and elucidated the underlying molecular mechanism. The present study opens the way for understanding the molecular mechanism of drug resistance through the direct quantitative comparison of ligand binding and enzymatic reaction with the same accuracy.

The Interplay of Electrostatics and Chemical Positioning in the Evolution of Antibiotic Resistance in TEM β-Lactamases

The interplay of enzyme active site electrostatics and chemical positioning is important for understanding the origin(s) of enzyme catalysis and the design of novel catalysts. We reconstruct the evolutionary trajectory of TEM-1 β-lactamase to TEM-52 toward extended-spectrum activity to better understand the emergence of antibiotic resistance and to provide insights into the structure-function paradigm and noncovalent interactions involved in catalysis. Utilizing a detailed kinetic analysis and the vibrational Stark effect, we quantify the changes in rates and electric fields in the Michaelis and acyl-enzyme complexes for penicillin G and cefotaxime to ascertain the evolutionary role of electric fields to modulate function. These data are combined with MD simulations to interpret and quantify the substrate-dependent structural changes during evolution. We observe that this evolutionary trajectory utilizes a large preorganized electric field and substrate-dependent chemical positioning to facilitate catalysis. This governs the evolvability, substrate promiscuity, and protein fitness landscape in TEM β-lactamase antibiotic resistance.

Free Energy Landscape and Proton Transfer Pathways of the Transimination Reaction at the Active site of the Serine Hydroxymethyltransferase Enzyme in Aqueous Medium

Serine hydroxymethyltransferase (SHMT) is a ubiquitous enzyme belonging to the fold type I or aspartate aminotransferase (AspAT) family of the pyridoxal 5'-phosphate (PLP)-dependent enzymes. Like other PLP-dependent enzymes, SHMT also undergoes the so-called transimination reaction before exhibiting its enzymatic activity. The transimination process constitutes an important pre-step for all PLP-dependent enzymes, where an internal aldimine of the PLP-enzyme complex gets converted to an external aldimine of the substrate-PLP complex at the active site of the enzyme. In case of the transimination reaction involving SHMT, the PLP molecule bound to the active site lysine residue of SHMT (internal aldimine) gets detached from the enzyme by a serine substrate to produce an external aldimine complex, where the PLP is now bound to the serine substrate. In the current study, the free energy surfaces and reaction pathways of different steps of the transimination reaction at the active site of SHMT are investigated by employing hybrid quantum mechanical/molecular mechanical (QM/MM) simulations combined with metadynamics methods of rare event sampling. It is found that the process of transimination involving serine and PLP at the active site of the SHMT enzyme takes place through different elementary steps such as the formation of the first geminal diamine intermediate (GDI1), transfer of a proton from the substrate serine to the phenolic oxygen of PLP, followed by another proton transfer from PLP to the amine nitrogen of lysine with the formation of the second geminal diamine intermediate (GDI2), and finally, detachment of the active site lysine residue from PLP to produce the external aldimine.

Probing the Suitability of Different Ca2+ Parameters for Long Simulations of Diisopropyl Fluorophosphatase

Organophosphate hydrolases are promising as potential biotherapeutic agents to treat poisoning with pesticides or nerve gases. However, these enzymes often need to be further engineered in order to become useful in practice. One example of such enhancement is the alteration of enantioselectivity of diisopropyl fluorophosphatase (DFPase). Molecular modeling techniques offer a unique opportunity to address this task rationally by providing a physical description of the substrate-binding process. However, DFPase is a metalloenzyme, and correct modeling of metal cations is a challenging task generally coming with a tradeoff between simulation speed and accuracy. Here, we probe several molecular mechanical parameter combinations for their ability to empower long simulations needed to achieve a quantitative description of substrate binding. We demonstrate that a combination of the Amber19sb force field with the recently developed 12-6 Ca2+ models allows us to both correctly model DFPase and obtain new insights into the DFP binding process.

Physics-based, neural network force fields for reactive molecular dynamics: Investigation of carbene formation from [EMIM+][OAc-]

Reactive molecular dynamics simulations enable a detailed understanding of solvent effects on chemical reaction mechanisms and reaction rates. While classical molecular dynamics using reactive force fields allows significantly longer simulation time scales and larger system sizes compared with ab initio molecular dynamics, constructing reactive force fields is a difficult and complex task. In this work, we describe a general approach following the empirical valence bond framework for constructing ab initio reactive force fields for condensed phase simulations by combining physics-based methods with neural networks (PB/NNs). The physics-based terms ensure the correct asymptotic behavior of electrostatic, polarization, and dispersion interactions and are compatible with existing solvent force fields. NNs are utilized for a versatile description of short-range orbital interactions within the transition state region and accurate rendering of vibrational motion of the reacting complex. We demonstrate our methodology for a simple deprotonation reaction of the 1-ethyl-3-methylimidazolium cation with acetate to form 1-ethyl-3-methylimidazol-2-ylidene and acetic acid. Our PB/NN force field exhibits ∼1 kJ mol^-1 mean absolute error accuracy within the transition state region for the gas-phase complex. To characterize the solvent modulation of the reaction profile, we compute potentials of mean force for the gas-phase reaction as well as the reaction within a four-ion cluster and benchmark against ab initio molecular dynamics simulations. We find that the surrounding ionic environment significantly destabilizes the formation of the carbene product, and we show that this effect is accurately captured by the reactive force field. By construction, the PB/NN potential may be directly employed for simulations of other solvents/chemical environments without additional parameterization.

On the Accuracy of QM/MM Models: A Systematic Study of Intramolecular Proton Transfer Reactions of Amino Acids in Water

This work presents a systematic assessment of QM/QM' and QM/MM models with respect to direct QM calculations for the tautomerization (neutral to zwitterion) reactions of amino acids (glycine, alanine, valine, aspartate, and neutral and protonated histidine) solvated in a 160 water cluster. The effect of varying QM region size and choice of embedding potentials, including fixed-charge and polarizable molecular mechanics force fields (TIP3P and EFP) and various semiempirical QM methods (PM7, GFN2-xTB, DFTBA, DFTB3, HF-3c, and PBEh-3c), on the accuracy of the models was examined. A surprising finding was that molecular mechanics force fields outperformed many of the semiempirical methods. Generally, the errors in the QM/QM' and QM/MM models converge slowly with respect to the QM region size, requiring 50 or more waters to be included in the QM region before the error in the model falls below 1 kcal mol^-1 of its pure QM result. Different QM region selection schemes were also compared, and it was found that selection based on Natural Population Analysis (NPA) atomic charges significantly reduced the error in the QM/QM' and QM/MM models particularly if a low-quality embedding potential was used. It is envisaged that these results will be useful for the development of future hybrid QM models.

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.

Cheminformatic quantum mechanical enzyme model design: A catechol-O-methyltransferase case study

To accurately simulate the inner workings of an enzyme active site with quantum mechanics (QM), not only must the reactive species be included in the model but also important surrounding residues, solvent, or coenzymes involved in crafting the microenvironment. Our lab has been developing the Residue Interaction Network Residue Selector (RINRUS) toolkit to utilize interatomic contact network information for automated, rational residue selection and QM-cluster model generation. Starting from an x-ray crystal structure of catechol-O-methyltransferase, RINRUS was used to construct a series of QM-cluster models. The reactant, product, and transition state of the methyl transfer reaction were computed for a total of 550 models, and the resulting free energies of activation and reaction were used to evaluate model convergence. RINRUS-designed models with only 200-300 atoms are shown to converge. RINRUS will serve as a cornerstone for improved and automated cheminformatics-based enzyme model design.

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.

Machine Learning for Chemical Reactions

Machine learning (ML) techniques applied to chemical reactions have a long history. The present contribution discusses applications ranging from small molecule reaction dynamics to computational platforms for reaction planning. ML-based techniques can be particularly relevant for problems involving both computation and experiments. For one, Bayesian inference is a powerful approach to develop models consistent with knowledge from experiments. Second, ML-based methods can also be used to handle problems that are formally intractable using conventional approaches, such as exhaustive characterization of state-to-state information in reactive collisions. Finally, the explicit simulation of reactive networks as they occur in combustion has become possible using machine-learned neural network potentials. This review provides an overview of the questions that can and have been addressed using machine learning techniques, and an outlook discusses challenges in this diverse and stimulating field. It is concluded that ML applied to chemistry problems as practiced and conceived today has the potential to transform the way with which the field approaches problems involving chemical reactions, in both research and academic teaching.

Comparative QM/MM Molecular Dynamics and Umbrella Sampling Simulations: Interaction of the Zinc-Bound Intermediate Gem-Diolate Trapoxin A Inhibitor and Acetyl-l-lysine Substrate with Histone Deacetylase 8

Targeting the genetic material without destruction is a priority to develop safe anticancer drugs. Histone deacetylase 8 (HDAC8), which is proved to be involved in carcinogenesis, is an enzyme associated with the chromatin for post-translational deacetylation of acetylated lysine. In this study, HDAC8 co-crystallized with the intermediate state tetrapeptide Trapoxin A (TA) inhibitor and the holoenzyme are utilized to find their conformational ensembles. Furthermore, the co-crystallized intermediate gem-diolate TA was used to find optimum interactions with the active site residues by conventional molecular dynamics (MD) simulation and QM/MM umbrella sampling. Finally, the intermediate state of the acetyl-l-lysine substrate was explored by QM/MM steered MD and compared to the binding of the intermediate state of the inhibitor. This research showed that HDAC8 is flexible and exists in conformational ensembles in its holoenzyme state. Binding of the intermediate state TA stabilizes its conformation. The optimum binding to the active site of HDAC8 for structures of gem-diolate TA (intermediate state) and acetyl-l-lysine (intermediate state) was determined according to the corresponding energy profiles. The use of these models will aid in the design of potentially reversible, potent, and selective inhibitors of HDAC8 for cancer treatment.

QM-Cluster Model Study of the Guaiacol Hydrogen Atom Transfer and Oxygen Rebound with Cytochrome P450 Enzyme GcoA

The key step of the O-demethylation of guaiacol by GcoA of the cytochrome P450-reductase pair was studied with DFT using two 10-residue and three 15-residue QM-cluster models. For each model, two reaction pathways were examined, beginning with a different guaiacol orientation. Based on this study, His354, Phe349, Glu249, and Pro250 residues were found to be important for keeping the heme in a planar geometry throughout the reaction. Val241 and Gly245 residues were needed in the QM-cluster models to provide the hydrophobic pocket for an appropriate guaiacol pose in the reaction. The aromatic triad Phe75, Phe169, and Phe395 may be necessary to facilitate guaiacol migrating into the enzyme active site, but it does not qualitatively affect kinetics and thermodynamics of the proposed mechanism. All QM-cluster models created by RINRUS agree very well with previous experimental work. This study provides details for better understanding enzymatic O-demethylation of lignins to form catechol derivatives by GcoA.

Biomolecular QM/MM Simulations: What Are Some of the "Burning Issues"?

QM/MM simulations have become an indispensable tool in many chemical and biochemical investigations. Considering the tremendous degree of success, including recognition by a 2013 Nobel Prize in Chemistry, are there still "burning challenges" in QM/MM methods, especially for biomolecular systems? In this short Perspective, we discuss several issues that we believe greatly impact the robustness and quantitative applicability of QM/MM simulations to many, if not all, biomolecules. We highlight these issues with observations and relevant advances from recent studies in our group and others in the field. Despite such limited scope, we hope the discussions are of general interest and will stimulate additional developments that help push the field forward in meaningful directions.

A QM/MM simulation study of transamination reaction at the active site of aspartate aminotransferase: Free energy landscape and proton transfer pathways

Transaminase is a key enzyme for amino acid metabolism, which reversibly catalyzes the transamination reaction with the help of PLP (pyridoxal 5^' -phosphate) as its cofactor. Here we have investigated the mechanism and free energy landscape of the transamination reaction involving the aspartate transaminase (AspTase) enzyme and aspartate-PLP (Asp-PLP) complex using QM/MM simulation and metadynamics methods. The reaction is found to follow a stepwise mechanism where the active site residue Lys258 acts as a base to shuttle a proton from α-carbon (CA) to imine carbon (C4A) of the PLP-Asp Schiff base. In the first step, the Lys258 abstracts the CA proton of the substrate leading to the formation of a carbanionic intermediate which is followed by the reprotonation of the Asp-PLP Schiff base at C4A atom by Lys258. It is found that the free energy barrier for the proton abstraction by Lys258 and that for the reprotonation are 17.85 and 3.57 kcal/mol, respectively. The carbanionic intermediate is 7.14 kcal/mol higher in energy than the reactant. Hence, the first step acts as the rate limiting step. The present calculations also show that the Lys258 residue undergoes a conformational change after the first step of transamination reaction and becomes proximal to C4A atom of the Asp-PLP Schiff base to favor the second step. The active site residues Tyr70* and Gly38 anchor the Lys258 in proper position and orientation during the first step of the reaction and stabilize the positive charge over Lys258 generated at the intermediate step.

Scalable molecular dynamics on CPU and GPU architectures with NAMD

NAMDis a molecular dynamics program designed for high-performance simulations of very large biological objects on CPU- and GPU-based architectures. NAMD offers scalable performance on petascale parallel supercomputers consisting of hundreds of thousands of cores, as well as on inexpensive commodity clusters commonly found in academic environments. It is written in C++ and leans on Charm++ parallel objects for optimal performance on low-latency architectures. NAMD is a versatile, multipurpose code that gathers state-of-the-art algorithms to carry out simulations in apt thermodynamic ensembles, using the widely popular CHARMM, AMBER, OPLS, and GROMOS biomolecular force fields. Here, we review the main features of NAMD that allow both equilibrium and enhanced-sampling molecular dynamics simulations with numerical efficiency. We describe the underlying concepts utilized by NAMD and their implementation, most notably for handling long-range electrostatics; controlling the temperature, pressure, and pH; applying external potentials on tailored grids; leveraging massively parallel resources in multiple-copy simulations; and hybrid quantum-mechanical/molecular-mechanical descriptions. We detail the variety of options offered by NAMD for enhanced-sampling simulations aimed at determining free-energy differences of either alchemical or geometrical transformations and outline their applicability to specific problems. Last, we discuss the roadmap for the development of NAMD and our current efforts toward achieving optimal performance on GPU-based architectures, for pushing back the limitations that have prevented biologically realistic billion-atom objects to be fruitfully simulated, and for making large-scale simulations less expensive and easier to set up, run, and analyze. NAMD is distributed free of charge with its source code at www.ks.uiuc.edu.

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.

A Review on Applications of Computational Methods in Drug Screening and Design

Drug development is one of the most significant processes in the pharmaceutical industry. Various computational methods have dramatically reduced the time and cost of drug discovery. In this review, we firstly discussed roles of multiscale biomolecular simulations in identifying drug binding sites on the target macromolecule and elucidating drug action mechanisms. Then, virtual screening methods (e.g., molecular docking, pharmacophore modeling, and QSAR) as well as structure- and ligand-based classical/de novo drug design were introduced and discussed. Last, we explored the development of machine learning methods and their applications in aforementioned computational methods to speed up the drug discovery process. Also, several application examples of combining various methods was discussed. A combination of different methods to jointly solve the tough problem at different scales and dimensions will be an inevitable trend in drug screening and design.

Automated exploration of the low-energy chemical space with fast quantum chemical methods

We propose and discuss an efficient scheme for the in silico sampling for parts of the molecular chemical space by semiempirical tight-binding methods combined with a meta-dynamics driven search algorithm. The focus of this work is set on the generation of proper thermodynamic ensembles at a quantum chemical level for conformers, but similar procedures for protonation states, tautomerism and non-covalent complex geometries are also discussed. The conformational ensembles consisting of all significantly populated minimum energy structures normally form the basis of further, mostly DFT computational work, such as the calculation of spectra or macroscopic properties. By using basic quantum chemical methods, electronic effects or possible bond breaking/formation are accounted for and a very reasonable initial energetic ranking of the candidate structures is obtained. Due to the huge computational speedup gained by the fast low-cost quantum chemical methods, overall short computation times even for systems with hundreds of atoms (typically drug-sized molecules) are achieved. Furthermore, specialized applications, such as sampling with implicit solvation models or constrained conformational sampling for transition-states, metal-, surface-, or noncovalently bound complexes are discussed, opening many possible applications in modern computational chemistry and drug discovery. The procedures have been implemented in a freely available computer code called CREST, that makes use of the fast and reliable GFNn-xTB methods.