Fast and flexible gpu accelerated binding free energy calculations within the amber molecular dynamics package

Alchemical free energy (AFE) calculations based on molecular dynamics (MD) simulations are key tools in both improving our understanding of a wide variety of biological processes and accelerating the design and optimization of therapeutics for numerous diseases. Computing power and theory have, however, long been insufficient to enable AFE calculations to be routinely applied in early stage drug discovery. One of the major difficulties in performing AFE calculations is the length of time required for calculations to converge to an ensemble average. CPU implementations of MD‐based free energy algorithms can effectively only reach tens of nanoseconds per day for systems on the order of 50,000 atoms, even running on massively parallel supercomputers. Therefore, converged free energy calculations on large numbers of potential lead compounds are often untenable, preventing researchers from gaining crucial insight into molecular recognition, potential druggability and other crucial areas of interest. Graphics Processing Units (GPUs) can help address this. We present here a seamless GPU implementation, within the PMEMD module of the AMBER molecular dynamics package, of thermodynamic integration (TI) capable of reaching speeds of >140 ns/day for a 44,907‐atom system, with accuracy equivalent to the existing CPU implementation in AMBER. The implementation described here is currently part of the AMBER 18 beta code and will be an integral part of the upcoming version 18 release of AMBER. © 2018 Wiley Periodicals, Inc.     

Theoretical study of morphine and heroin: Conformational study in gas phase and aqueous solution and electron distribution analysis

The conformational preferences of morphine and heroin were studied in gas phase and with inclusion of solvent effects. At 298.15 K, three conformers are significant for isolated morphine, all of them displaying antiperiplanar arrangement for the C2? C3? O? H unit, and there is only one significantly populated conformer for heroin. Quantum theory of atoms in molecules analysis of the electron density in their most populated conformers in gas phase indicates that the positive charge is shared among the amino hydrogen, those hydrogens of the methylamino group, and all of the hydrogens attached to the bridgehead carbons. © 2010 Wiley Periodicals, Inc. Int J Quantum Chem 110:2472–2482, 2010     

Electron density analysis of large (molecular and periodic) systems: A parallel implementation

A parallel implementation is presented of a series of algorithms for the evaluation of several one‐electron properties of large molecular and periodic (of any dimensionality) systems. The electron charge and momentum densities of the system, the electrostatic potential, X‐ray structure factors, directional Compton profiles can be effectively evaluated at low computational cost along with a full topological analysis of the electron charge density (ECD) of the system according to Bader's quantum theory of atoms in molecules. The speedup of the parallelization of the different algorithms is presented. The search of all symmetry‐irreducible critical points of the ECD of the crystallized crambin protein and the evaluation of all the corresponding bond paths, for instance, would require about 32 days if run in serial mode and reduces to less than 2 days when run in parallel mode over 32 processors. © 2015 Wiley Periodicals, Inc.     

A new set of atomic radii for accurate estimation of solvation free energy by Poisson–Boltzmann solvent model

The Poisson–Boltzmann implicit solvent (PB) is widely used to estimate the solvation free energies of biomolecules in molecular simulations. An optimized set of atomic radii (PB radii) is an important parameter for PB calculations, which determines the distribution of dielectric constants around the solute. We here present new PB radii for the AMBER protein force field to accurately reproduce the solvation free energies obtained from explicit solvent simulations. The presented PB radii were optimized using results from explicit solvent simulations of the large systems. In addition, we discriminated PB radii for N‐ and C‐terminal residues from those for nonterminal residues. The performances using our PB radii showed high accuracy for the estimation of solvation free energies at the level of the molecular fragment. The obtained PB radii are effective for the detailed analysis of the solvation effects of biomolecules. © 2014 The Authors Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

An ant algorithm for the conformational analysis of flexible molecules

Originally, the ant system was developed for optimization in discrete search spaces such as the traveling salesman problem. We detail our adaptation of the algorithm to optimization in the continuous search space of conformational analysis. The parameters of the algorithm were tuned using a simple test molecule, undecane, and a drug molecule, imatinib. The algorithm is further tested on four more drug or drug-like molecules, on vitamin A and on alanine tetrapeptide.     

Localization and quantification of hydrophobicity: The molecular free energy density (MolFESD) concept and its application to sweetness recognition

A method for the localization, the quantification, and the analysis of hydrophobicity of a molecule or a molecular fragment is presented. It is shown that the free energy of solvation for a molecule or the transfer free energy from one solvent to another can be represented by a surface integral of a scalar quantity, the molecular free energy surface density (MolFESD), over the solvent accessible surface of that molecule. This MolFESD concept is based on a model approach where the solvent molecules are considered to be small in comparison to the solute molecule, and the solvent can be represented by a continuous medium with a given dielectric constant. The transfer energy surface density for a 1-octanol/water system is empirically determined employing a set of atomic increment contributions and distance dependent membership functions measuring the contribution of the increments to the surface value of the MolFESD. The MolFESD concept can be well used for the quantification of the purely hydrophobic contribution to the binding constants of molecule-receptor complexes. This is demonstrated with the sweeteners sucrose and sucralose and various halogen derivatives. Therein the relative sweetness, which is assumed to be proportional to the binding constant, nicely correlates to the surface integral over the positive, hydrophobic part of the MolFESD, indicating that the sweetness receptor can be characterized by a highly flexible hydrophobic pocket instead of a localized binding site.     

ERKALE—A flexible program package for X‐ray properties of atoms and molecules

ERKALE is a novel software program for computing X‐ray properties, such as ground‐state electron momentum densities, Compton profiles, and core and valence electron excitation spectra of atoms and molecules. The program operates at Hartree–Fock or density‐functional level of theory and supports Gaussian basis sets of arbitrary angular momentum and a wide variety of exchange‐correlation functionals. ERKALE includes modern convergence accelerators such as Broyden and ADIIS and it is suitable for general use, as calculations with thousands of basis functions can routinely be performed on desktop computers. Furthermore, ERKALE is written in an object oriented manner, making the code easy to understand and to extend to new properties while being ideal also for teaching purposes. © 2012 Wiley Periodicals, Inc.     

BAND NN: A Deep Learning Framework for Energy Prediction and Geometry Optimization of Organic Small Molecules

Recent advances in artificial intelligence along with the development of large data sets of energies calculated using quantum mechanical (QM)/density functional theory (DFT) methods have enabled prediction of accurate molecular energies at reasonably low computational cost. However, machine learning models that have been reported so far require the atomic positions obtained from geometry optimizations using high-level QM/DFT methods as input in order to predict the energies and do not allow for geometry optimization. In this study, a transferable and molecule size-independent machine learning model bonds (B), angles (A), nonbonded (N) interactions, and dihedrals (D) neural network (BAND NN) based on a chemically intuitive representation inspired by molecular mechanics force fields is presented. The model predicts the atomization energies of equilibrium and nonequilibrium structures as sum of energy contributions from bonds (B), angles (A), nonbonds (N), and dihedrals (D) at remarkable accuracy. The robustness of the proposed model is further validated by calculations that span over the conformational, configurational, and reaction space. The transferability of this model on systems larger than the ones in the data set is demonstrated by performing calculations on selected large molecules. Importantly, employing the BAND NN model, it is possible to perform geometry optimizations starting from nonequilibrium structures along with predicting their energies. © 2019 Wiley Periodicals, Inc.     

Simulating GTP:Mg and GDP:Mg with a simple force field: A structural and thermodynamic analysis

Di‐ and tri‐phosphate nucleotides are essential cofactors for many proteins, usually in an Mg²⁺‐bound form. Proteins like GTPases often detect the difference between NDP and NTP and respond by changing conformations. To study such complexes, simple, fixed charge force fields have been used, which allow long simulations and precise free energy calculations. The preference for NTP or NDP binding depends on many factors, including ligand structure and Mg²⁺ coordination and the changes they undergo upon binding. Here, we use a simple force field to examine two Mg²⁺ coordination modes for the unbound GDP and GTP: direct, or “Inner Sphere” (IS) coordination by one or more phosphate oxygens and indirect, “Outer Sphere” (OS) coordination involving one or more bridging waters. We compare GTP: and GDP:Mg binding with OS and IS coordination; combining the results with experimental data then indicates that GTP prefers the latter. We also examine different kinds of IS coordination and their sensitivity to a key force field parameter: the optimal Mg:oxygen van der Waals distance R_min. Increasing R_min improves the Mg:oxygen distances, the GTP: and GDP:Mg binding affinities, and the fraction of GTP:Mg with β + γ phosphate coordination, but does not improve or change the GTP/GDP affinity difference, which remains much larger than experiment. It has no effect on the free energy of GDP binding to a GTPase. © 2012 Wiley Periodicals, Inc.     

A new formulation of the PCM solvation method: PCM-QINTn

A new formulation of the PCM electrostatic solution problem is proposed. Through a new derivation of the PCM-CLSn expression we propose an interpolation formula that improved the convergence: PCM-QINTn. All the available formulations are applied to the evaluation of the electrostatic component of the free energy of solvation for some molecular systems. In addition, PCM-QINT derivatives of G _el with respect to atomic coordinates are evaluated. The computational costs are compared with those of PCM-direct formulation. Received: 21 October 1996 / Accepted: 7 January 1997     

Heuristic molecular lipophilicity potential (HMLP): lipophilicity and hydrophilicity of amino acid side chains

Heuristic molecular lipophilicity potential (HMLP) is applied in the study of lipophilicity and hydrophilicity of 20 natural amino acids side chains. The HMLP parameters, surface area S(i), lipophilic indices L(i), and hydrophilic indices H(i) of amino acid side chains are derived from lipophilicity potential L(r). The parameters are correlated with the experimental data of phase-transferring free energies of vapor-to-water, vapor-to-cyclohexane, vapor-to-octanol, cyclohexane-to-water, octanol-to-water, and cyclohexane-to-octanol through a linear free energy equation DeltaG(0)(tr,i) = b(0) + b(1)S(i) (+) + b(2)S(i) (-) + b(3)L(i) + b(4)H(i). For all above six phase-transfer free energies, the HMLP parameters of 20 amino acid side chains give good calculation results using linear free energy equation. HMLP is an ab initio quantum chemical approach and a structure-based technique. Except for atomic van der Waals radii, there are no other empirical parameters used. The HMLP has clear physical and chemical meaning and provides useful lipophilic and hydrophilic parameters for the studies of proteins and peptides.     

A conformational analysis method for understanding the energy landscapes of clusters.

A newly developed unbiased structural optimization method, named dynamic lattice searching (DLS), is proposed as an approach for conformational analysis of atomic/molecular clusters and used in understanding the energy landscape of large clusters. The structures of clusters are described in terms of the number of basic tetrahedron (BT) units they contain. We found that the hit numbers of different structural motifs in DLS runs is proportional to the number of BTs. A parameter T(max) is defined to limit the maximal number of atoms moved in a structural transition. Results show that T(max) is a key parameter for modulating the efficiency of the DLS method and has a great influence on the hit number of different motifs in DLS runs. Finally, the effect of potential range on the conformational distribution of the (Morse)(98) cluster is also discussed with different potential-range parameters.     

Toward quantitative estimates of binding affinities for protein–ligand systems involving large inhibitor compounds: A steered molecular dynamics simulation route

Understanding binding mechanisms between enzymes and potential inhibitors and quantifying protein – ligand affinities in terms of binding free energy is of primary importance in drug design studies. In this respect, several approaches based on molecular dynamics simulations, often combined with docking techniques, have been exploited to investigate the physicochemical properties of complexes of pharmaceutical interest. Even if the geometric properties of a modeled protein – ligand complex can be well predicted by computational methods, it is still challenging to rank with chemical accuracy a series of ligand analogues in a consistent way. In this article, we face this issue calculating relative binding free energies of a focal adhesion kinase, an important target for the development of anticancer drugs, with pyrrolopyrimidine‐based ligands having different inhibitory power. To this aim, we employ steered molecular dynamics simulations combined with nonequilibrium work theorems for free energy calculations. This technique proves very powerful when a series of ligand analogues is considered, allowing one to tackle estimation of protein – ligand relative binding free energies in a reasonable time. In our cases, the calculated binding affinities are comparable with those recovered from experiments by exploiting the Michaelis – Menten mechanism with a competitive inhibitor.     

A novel approach to the calculation of the free energy due to molecular flexibility

The contribution of the molecular flexibility to the solvation excess free energy is expressed in terms of probabilities of reaching hard limits on intramolecular coordinates in a series of calculations successively relaxing those limits. Numerical tests on the harmonic oscillator are also presented and used to make suggestion about computational issues.     

A molecular mechanics model for imatinib and imatinib:kinase binding

Imatinib is an important anticancer drug, which binds specifically to the Abl kinase and blocks its signalling activity. To model imatinib:protein interactions, we have developed a molecular mechanics force field for imatinib and four close analogues, which is consistent with the CHARMM force field for proteins and nucleic acids. Atomic charges and Lennard‐Jones parameters were derived from a supermolecule ab initio approach. We considered the ab initio energies and geometries of a probe water molecule interacting with imatinib fragments at 32 different positions. We considered both a neutral and a protonated imatinib. The final RMS deviation between the ab initio and force field energies, averaged over both forms, was 0.2 kcal/mol. The model also reproduces the ab initio geometry and flexibility of imatinib. To apply the force field to imatinib:Abl simulations, it is also necessary to determine the most likely imatinib protonation state when it binds to Abl. This was done using molecular dynamics free energy simulations, where imatinib is reversibly protonated during a series of MD simulations, both in solution and in complex with Abl. The simulations indicate that imatinib binds to Abl in its protonated, positively‐charged form. To help test the force field and the protonation prediction, we did MD free energy simulations that compare the Abl binding affinities of two imatinib analogs, obtaining good agreement with experiment. Finally, two new imatinib variants were considered, one of which is predicted to have improved Abl binding. This variant could be of interest as a potential drug. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Free energies of solvation in the context of protein folding: Implications for implicit and explicit solvent models

Implicit solvent models for biomolecular simulations have been developed to use in place of more expensive explicit models; however, these models make many assumptions and approximations that are likely to affect accuracy. Here, the changes in free energies of solvation upon folding of several fast folding proteins are calculated from previously run μs–ms simulations with a number of implicit solvent models and compared to the values needed to be consistent with the explicit solvent model used in the simulations. In the majority of cases, there is a significant and substantial difference between the values calculated from the two approaches that is robust to the details of the calculations. These differences could only be remedied by selecting values for the model parameters—the internal dielectric constant for the polar term and the surface tension coefficient for the nonpolar term—that were system‐specific or physically unrealistic. We discuss the potential implications of our findings for both implicit and explicit solvent simulations. © 2015 Wiley Periodicals, Inc.