RESUMO
We present a comprehensive study investigating the potential gain in accuracy for calculating absolute solvation free energies (ASFE) using a neural network potential to describe the intramolecular energy of the solute. We calculated the ASFE for most compounds from the FreeSolv database using the Open Force Field (OpenFF) and compared them to earlier results obtained with the CHARMM General Force Field (CGenFF). By applying a nonequilibrium (NEQ) switching approach between the molecular mechanics (MM) description (either OpenFF or CGenFF) and the neural net potential (NNP)/MM level of theory (using ANI-2x as the NNP potential), we attempted to improve the accuracy of the calculated ASFEs. The predictive performance of the results did not change when this approach was applied to all 589 small molecules in the FreeSolv database that ANI-2x can describe. When selecting a subset of 156 molecules, focusing on compounds where the force fields performed poorly, we saw a slight improvement in the root-mean-square error (RMSE) and mean absolute error (MAE). The majority of our calculations utilized unidirectional NEQ protocols based on Jarzynski's equation. Additionally, we conducted bidirectional NEQ switching for a subset of 156 solutes. Notably, only a small fraction (10 out of 156) exhibited statistically significant discrepancies between unidirectional and bidirectional NEQ switching free energy estimates.
RESUMO
Double decoupling absolute binding free energy simulations require an intermediate state at which the ligand is held solely by restraints in a position and orientation resembling the bound state. One possible choice consists of one distance, two angle, and three dihedral angle restraints. Here, I demonstrate that in practically all cases the analytical correction derived under the rigid rotator harmonic oscillator approximation is sufficient to account for the free energy of the restraints.
Assuntos
Termodinâmica , Ligantes , Ligação Proteica , Proteínas/química , Proteínas/metabolismo , Simulação de Dinâmica Molecular , Modelos MolecularesRESUMO
To achieve chemical accuracy in free energy calculations, it is necessary to accurately describe the system's potential energy surface and efficiently sample configurations from its Boltzmann distribution. While neural network potentials (NNPs) have shown significantly higher accuracy than classical molecular mechanics (MM) force fields, they have a limited range of applicability and are considerably slower than MM potentials, often by orders of magnitude. To address this challenge, Rufa et al. [Rufa et al. bioRxiv 2020, 10.1101/2020.07.29.227959.] suggested a two-stage approach that uses a fast and established MM alchemical energy protocol, followed by reweighting the results using NNPs, known as endstate correction or indirect free energy calculation. This study systematically investigates the accuracy and robustness of reweighting from an MM reference to a neural network target potential (ANI-2x) for an established data set in vacuum, using single-step free-energy perturbation (FEP) and nonequilibrium (NEQ) switching simulation. We assess the influence of longer switching lengths and the impact of slow degrees of freedom on outliers in the work distribution and compare the results to those of multistate equilibrium free energy simulations. Our results demonstrate that free energy calculations between NNPs and MM potentials should be preferably performed using NEQ switching simulations to obtain accurate free energy estimates. NEQ switching simulations between the MM potentials and NNPs are efficient, robust, and trivial to implement.
RESUMO
We recently introduced transformato, an open-source Python package for the automated setup of large-scale calculations of relative solvation and binding free energy differences. Here, we extend the capabilities of transformato to the calculation of absolute solvation free energy differences. After careful validation against the literature results and reference calculations with the PERT module of CHARMM, we used transformato to compute absolute solvation free energies for most molecules in the FreeSolv database (621 out of 642). The force field parameters were obtained with the program cgenff (v2.5.1), which derives missing parameters from the CHARMM general force field (CGenFF v4.6). A long-range correction for the Lennard-Jones interactions was added to all computed solvation free energies. The mean absolute error compared to the experimental data is 1.12 kcal/mol. Our results allow a detailed comparison between the AMBER and CHARMM general force fields and provide a more in-depth understanding of the capabilities and limitations of the CGenFF small molecule parameters.
RESUMO
Non-equilibrium work switching simulations and Jarzynski's equation are a reliable method for computing free energy differences, ΔAlowâhigh, between two levels of theory, such as a pure molecular mechanical (MM) and a quantum mechanical/molecular mechanical (QM/MM) description of a system of interest. Despite the inherent parallelism, the computational cost of this approach can quickly become very high. This is particularly true for systems where the core region, the part of the system to be described at different levels of theory, is embedded in an environment such as explicit solvent water. We find that even for relatively simple solute-water systems, switching lengths of at least 5 ps are necessary to compute ΔAlowâhigh reliably. In this study, we investigate two approaches towards an affordable protocol, with an emphasis on keeping the switching length well below 5 ps. Inserting a hybrid charge intermediate state with modified partial charges, which resembles the charge distribution of the desired high level, makes it possible to obtain reliable calculations with 2 ps switches. Attempts using step-wise linear switching paths, on the other hand, did not lead to improvement, i.e., a faster convergence for all systems. To understand these findings, we analyzed the solutes' properties as a function of the partial charges used and the number of water molecules in direct contact with the solute, and studied the time needed for water molecules to reorient themselves upon a change in the solute's charge distribution.
RESUMO
We present the software package transformato for the setup of large-scale relative binding free energy calculations. Transformato is written in Python as an open source project (https://github.com/wiederm/transformato); in contrast to comparable tools, it is not closely tied to a particular molecular dynamics engine to carry out the underlying simulations. Instead of alchemically transforming a ligand L 1 directly into another L 2, the two ligands are mutated to a common core. Thus, while dummy atoms are required at intermediate states, in particular at the common core state, none are present at the physical endstates. To validate the method, we calculated 76 relative binding free energy differences Δ Δ G L 1 â L 2 b i n d for five protein-ligand systems. The overall root mean squared error to experimental binding free energies is 1.17 kcal/mol with a Pearson correlation coefficient of 0.73. For selected cases, we checked that the relative binding free energy differences between pairs of ligands do not depend on the choice of the intermediate common core structure. Additionally, we report results with and without hydrogen mass reweighting. The code currently supports OpenMM, CHARMM, and CHARMM/OpenMM directly. Since the program logic to choose and construct alchemical transformation paths is separated from the generation of input and topology/parameter files, extending transformato to support additional molecular dynamics engines is straightforward.
RESUMO
We describe the theory of the so-called common-core/serial-atom-insertion (CC/SAI) approach to compute alchemical free energy differences and its practical implementation in a Python package called Transformato. CC/SAI is not tied to a specific biomolecular simulation program and does not rely on special purpose code for alchemical transformations. To calculate the alchemical free energy difference between several small molecules, the physical end-states are mutated into a suitable common core. Since this only requires turning off interactions, the setup of intermediate states is straightforward to automate. Transformato currently supports CHARMM and OpenMM as back ends to carry out the necessary molecular dynamics simulations, as well as post-processing calculations. We validate the method by computing a series of relative solvation free energy differences.
Assuntos
Simulação de Dinâmica Molecular , Entropia , TermodinâmicaRESUMO
A key step during indirect alchemical free energy simulations using quantum mechanical/molecular mechanical (QM/MM) hybrid potential energy functions is the calculation of the free energy difference ΔAlowâhigh between the low level (e.g., pure MM) and the high level of theory (QM/MM). A reliable approach uses nonequilibrium work (NEW) switching simulations in combination with Jarzynski's equation; however, it is computationally expensive. In this study, we investigate whether it is more efficient to use more shorter switches or fewer but longer switches. We compare results obtained with various protocols to reference free energy differences calculated with Crooks' equation. The central finding is that fewer longer switches give better converged results. As few as 200 sufficiently long switches lead to ΔAlowâhigh values in good agreement with the reference results. This optimized protocol reduces the computational cost by a factor of 40 compared to earlier work. We also describe two tools/ways of analyzing the raw data to detect sources of poor convergence. Specifically, we find it helpful to analyze the raw data (work values from the NEW switching simulations) in a quasi-time series-like manner. Principal component analysis helps to detect cases where one or more conformational degrees of freedom are different at the low and high level of theory.
Assuntos
Teoria Quântica , Entropia , Conformação Molecular , TermodinâmicaRESUMO
In calculations of relative free energy differences, the number of atoms of the initial and final states is rarely the same. This necessitates the introduction of dummy atoms. These placeholders interact with the physical system only by bonded energy terms. We investigate the conditions necessary so that the presence of dummy atoms does not influence the result of a relative free energy calculation. On the one hand, one has to ensure that dummy atoms only give a multiplicative contribution to the partition function so that their contribution cancels from double-free energy differences. On the other hand, the bonded terms used to attach a dummy atom (or group of dummy atoms) to the physical system have to maintain it in a well-defined position and orientation relative to the physical system. A detailed theoretical analysis of both aspects is provided, illustrated by 24 calculations of relative solvation free energy differences, for which all four legs of the underlying thermodynamic cycle were computed. Cycle closure (or lack thereof) was used as a sensitive indicator to probing the effects of dummy atom treatment on the resulting free energy differences. We find that a naive (but often practiced) treatment of dummy atoms results in errors of up to kBT when calculating the relative solvation free energy difference between two small solutes, such as methane and ammonia. While our analysis focuses on the so-called single topology approach to set up alchemical transformations, similar considerations apply to dual topology, at least many widely used variants thereof.
RESUMO
Ionic liquids are an interesting class of soft matter with viscosities of one or two orders of magnitude higher than that of water. Unfortunately, classical, non-polarizable molecular dynamics (MD) simulations of ionic liquids result in too slow dynamics and demonstrate the need for explicit inclusion of polarizability. The inclusion of polarizability, here via the Drude oscillator model, requires amendments to the employed thermostat, where we consider a dual Nosé-Hoover thermostat, as well as a dual Langevin thermostat. We investigate the effects of the choice of a thermostat and the underlying parameters such as the masses and force constants of the Drude particles on static and dynamic properties of ionic liquids. Here, we show that Langevin thermostats are not suitable for investigating the dynamics of ionic liquids. Since polarizable MD simulations are associated with high computational costs, we employed a self-developed graphics processing unit enhanced code within the MD program CHARMM to keep the overall computational effort reasonable.
RESUMO
The use of the most accurate (i.e., QM or QM/MM) levels of theory for free energy simulations (FES) is typically not possible. Primarily, this is because the computational cost associated with the extensive configurational sampling needed for converging FES is prohibitive. To ensure the feasibility of QM-based FES, the "indirect" approach is generally taken, necessitating a free energy calculation between the MM and QM/MM potential energy surfaces. Ideally, this step is performed with standard free energy perturbation (Zwanzig's equation) as it only requires simulations be carried out at the low level of theory; however, work from several groups over the past few years has conclusively shown that Zwanzig's equation is ill-suited to this task. As such, many approximations have arisen to mitigate difficulties with Zwanzig's equation. One particularly popular notion is that the convergence of Zwanzig's equation can be improved by using interaction energy differences instead of total energy differences. Although problematic numerical fluctuations (a major problem when using Zwanzig's equation) are indeed reduced, our results and analysis demonstrate that this "interaction energy approximation" (IEA) is theoretically incorrect, and the implicit approximation invoked is spurious at best. Herein, we demonstrate this via solvation free energy calculations using IEA from two different low levels of theory to the same target high level. Results from this proof-of-concept consistently yield the wrong results, deviating by â¼1.5 kcal/mol from the rigorously obtained value.
RESUMO
Indirect (S)QM/MM free energy simulations (FES) are vital to efficiently incorporating sufficient sampling and accurate (QM) energetic evaluations when estimating free energies of practical/experimental interest. Connecting between levels of theory, i.e., calculating Δ A l o w â h i g h , remains to be the most challenging step within an indirect FES protocol. To improve calculations of Δ A l o w â h i g h , we must: (1) compare the performance of all FES methods currently available; and (2) compile and maintain datasets of Δ A l o w â h i g h calculated for a wide-variety of molecules so that future practitioners may replicate or improve upon the current state-of-the-art. Towards these two aims, we introduce a new dataset, "HiPen", which tabulates Δ A g a s M M â 3 o b (the free energy associated with switching from an M M to an S C C - D F T B molecular description using the 3ob parameter set in gas phase), calculated for 22 drug-like small molecules. We compare the calculation of this value using free energy perturbation, Bennett's acceptance ratio, Jarzynski's equation, and Crooks' equation. We also predict the reliability of each calculated Δ A g a s M M â 3 o b by evaluating several convergence criteria including sample size hysteresis, overlap statistics, and bias metric ( Π ). Within the total dataset, three distinct categories of molecules emerge: the "good" molecules, for which we can obtain converged Δ A g a s M M â 3 o b using Jarzynski's equation; "bad" molecules which require Crooks' equation to obtain a converged Δ A g a s M M â 3 o b ; and "ugly" molecules for which we cannot obtain reliably converged Δ A g a s M M â 3 o b with either Jarzynski's or Crooks' equations. We discuss, in depth, results from several example molecules in each of these categories and describe how dihedral discrepancies between levels of theory cause convergence failures even for these gas phase free energy simulations.
Assuntos
Metabolismo Energético , Proteínas/metabolismo , Termodinâmica , Água/metabolismo , Entropia , Teoria QuânticaRESUMO
The calculation of free energy differences between levels of theory has numerous potential pitfalls. Chief among them is the lack of overlap, i.e., ensembles generated at one level of theory (e.g., "low") not being good approximations of ensembles at the other (e.g., "high"). Numerous strategies have been devised to mitigate this issue. However, the most straightforward approach is to ensure that the "low" level ensemble more closely resembles that of the "high". Ideally, this is done without increasing computational cost. Herein, we demonstrate that by reparametrizing classical intramolecular potentials to reproduce high level forces (i.e., force matching) configurational overlap between a "low" (i.e., classical) and "high" (i.e., quantum) level can be significantly improved. This procedure is validated on two test cases and results in vastly improved convergence of free energy simulations.
RESUMO
The structural resolution of a bound ligand-receptor complex is a key asset to efficiently drive lead optimization in drug design. However, structural resolution of many drug targets still remains a challenging endeavor. In the absence of structural knowledge, scientists resort to structure-activity relationships (SARs) to promote compound development. In this study, we incorporated ligand-based knowledge to formulate a docking scoring function that evaluates binding poses for their agreement with a known SAR. We showcased this protocol by identifying the binding mode of the pyrazoloquinolinone (PQ) CGS-8216 at the benzodiazepine binding site of the GABAA receptor. Further evaluation of the final pose by molecular dynamics and free energy simulations revealed a close proximity between the pendent phenyl ring of the PQ and γ2D56, congruent with the low potency of carboxyphenyl analogues. Ultimately, we introduced the γ2D56A mutation and in fact observed a 10-fold potency increase in the carboxyphenyl analogue, providing experimental evidence in favor of our binding hypothesis.
Assuntos
Pirazóis/farmacologia , Receptores de GABA-A/metabolismo , Benzodiazepinas/metabolismo , Sítios de Ligação , Humanos , Ligantes , Simulação de Acoplamento Molecular , Subunidades Proteicas/química , Subunidades Proteicas/metabolismo , Pirazóis/química , Receptores de GABA-A/química , Software , Relação Estrutura-AtividadeRESUMO
We demonstrate that Jarzynski's equation can be used to reliably compute free energy differences between low and high level representations of systems. The need for such a calculation arises when employing the so-called "indirect" approach to free energy simulations with mixed quantum mechanical/molecular mechanical (QM/MM) Hamiltonians; a popular technique for circumventing extensive simulations involving quantum chemical computations. We have applied this methodology to several small and medium sized organic molecules, both in the gas phase and explicit solvent. Test cases include several systems for which the standard approach; that is, free energy perturbation between low and high level description, fails to converge. Finally, we identify three major areas in which the difference between low and high level representations make the calculation of ΔAlowâhigh difficult: bond stretching and angle bending, different preferred conformations, and the response of the MM region to the charge distribution of the QM region. © 2016 Wiley Periodicals, Inc.
RESUMO
We present a new approach that incorporates flexibility based on extensive MD simulations of protein-ligand complexes into structure-based pharmacophore modeling and virtual screening. The approach uses the multiple coordinate sets saved during the MD simulations and generates for each frame a pharmacophore model. Pharmacophore models with the same pharmacophore features are pooled. In this way the high number of pharmacophore models that results from the MD simulation is reduced to only a few hundred representative pharmacophore models. Virtual screening runs are performed with every representative pharmacophore model; the screening results are combined and rescored to generate a single hit-list. The score for a particular molecule is calculated based on the number of representative pharmacophore models which classified it as active. Hence, the method is called common hits approach (CHA). The steps between the MD simulation and the final hit-list are performed automatically and without user interaction. We test the performance of CHA for virtual screening using screening databases with active and inactive compounds for 40 protein-ligand systems. The results of the CHA are compared to the (i) median screening performance of all representative pharmacophore models of protein-ligand systems, as well as to the virtual screening performance of (ii) a random classifier, (iii) the pharmacophore model derived from the experimental structure in the PDB, and (iv) the representative pharmacophore model appearing most frequently during the MD simulation. For the 34 (out of 40) protein-ligand complexes, for which at least one of the approaches was able to perform better than a random classifier, the highest enrichment was achieved using CHA in 68% of the cases, compared to 12% for the PDB pharmacophore model and 20% for the representative pharmacophore model appearing most frequently. The availabilithy of diverse sets of different pharmacophore models is utilized to analyze some additional questions of interest in 3D pharmacophore-based virtual screening.
Assuntos
Avaliação Pré-Clínica de Medicamentos/métodos , Simulação de Dinâmica Molecular , Ligantes , Proteínas/química , Proteínas/metabolismo , Interface Usuário-ComputadorRESUMO
ABSTRACT: Pharmacophore modeling is a widely used technique in computer-aided drug discovery. Structure-based pharmacophore models of a ligand in complex with a protein have proven to be useful for supporting in silico hit discovery, hit to lead expansion, and lead optimization. As a structure-based approach it depends on the correct interpretation of ligand-protein interactions. There are legitimate concerns about the fidelity of the bound ligand and about non-physiological contacts with parts of the crystal and the solvent effects that influence the protein structure. A possible way to refine the structure of a protein-ligand system is to use the final structure of a given MD simulation. In this study we compare pharmacophore models built using the initial protein-ligand structure obtained from the protein data bank (PDB) with pharmacophore models built with the final structure of a molecular dynamics simulation. We show that the pharmacophore models differ in feature number and feature type and that the pharmacophore models built from the last structure of a MD simulation shows in some cases better ability to distinguish between active and decoy ligand structures.
RESUMO
Molecular dynamics simulations of twelve protein-ligand systems were used to derive a single, structure based pharmacophore model for each system. These merged models combine the information from the initial experimental structure and from all snapshots saved during the simulation. We compared the merged pharmacophore models with the corresponding PDB pharmacophore models, i.e., the static models generated from an experimental structure in the usual manner. The frequency of individual features, of feature types and the occurrence of features not present in the static model derived from the experimental structure were analyzed. We observed both pharmacophore features not visible in the traditional approach, as well as features which disappeared rapidly during the molecular dynamics simulations and which may well be artifacts of the initial PDB structure-derived pharmacophore model. Our approach helps mitigate the sensitivity of structure based pharmacophore models to the single set of coordinates present in the experimental structure. Further, the frequency with which specific features occur during the MD simulation may aid in ranking the importance of individual features.
Assuntos
Desenho de Fármacos , Modelos Químicos , Simulação de Dinâmica Molecular , Proteínas/química , Proteínas/ultraestrutura , Sítios de Ligação , Ligação de Hidrogênio , Ligantes , Ligação Proteica , Conformação ProteicaRESUMO
Carrying out free energy simulations (FES) using quantum mechanical (QM) Hamiltonians remains an attractive, albeit elusive goal. Renewed efforts in this area have focused on using "indirect" thermodynamic cycles to connect "low level" simulation results to "high level" free energies. The main obstacle to computing converged free energy results between molecular mechanical (MM) and QM (ΔA(MMâQM)), as recently demonstrated by us and others, is differences in the so-called "stiff" degrees of freedom (e.g., bond stretching) between the respective energy surfaces. Herein, we demonstrate that this problem can be efficiently circumvented using nonequilibrium work (NEW) techniques, i.e., Jarzynski's and Crooks' equations. Initial applications of computing ΔA(NEW)(MMâQM), for blocked amino acids alanine and serine as well as to generate butane's potentials of mean force via the indirect QM/MM FES method, showed marked improvement over traditional FES approaches.
RESUMO
BACKGROUND: Accurately modeling condensed phase processes is one of computation's most difficult challenges. Include the possibility that conformational dynamics may be coupled to chemical reactions, where multiscale (i.e., QM/MM) methods are needed, and this task becomes even more daunting. METHODS: Free energy simulations (i.e., molecular dynamics), multiscale modeling, and reweighting schemes. RESULTS: Herein, we present two new approaches for mitigating the aforementioned challenges. The first is a new chain-of-replica method (off-path simulations, OPS) for computing potentials of mean force (PMFs) along an easily defined reaction coordinate. This development is coupled with a new distributed, highly-parallel replica framework (REPDstr) within the CHARMM package. Validation of these new schemes is carried out on two processes that undergo conformational changes. First is the simple torsional rotation of butane, while a much more challenging glycosidic rotation (in vacuo and solvated) is the second. Additionally, a new approach that greatly improves (i.e., possibly an order of magnitude) the efficiency of computing QM/MM PMFs is introduced and compared to standard schemes. Our efforts are grounded in the recently developed method for efficiently computing QM-based free energies (i.e., QM-Non-Boltzmann Bennett, QM-NBB). Again, we validate this new technique by computing the QM/MM PMF of butane's torsional rotation. CONCLUSIONS: The OPS-REPDstr method is a promising new approach that overcomes many limitations of standard pathway simulations in CHARMM. The combination of QM-NBB with pathway techniques is very promising as it offers significant advantages over current procedures. GENERAL SIGNIFICANCE: Efficiently computing potentials of mean force is a major, unresolved, area of interest. This article is part of a Special Issue entitled Recent developments of molecular dynamics.