Free energy perturbation calculations involving potential function changes

Alchemical Binding Free Energy Calculations in AMBER20: Advances and Best Practices for Drug Discovery

Predicting protein-ligand binding affinities and the associated thermodynamics of biomolecular recognition is a primary objective of structure-based drug design. Alchemical free energy simulations offer a highly accurate and computationally efficient route to achieving this goal. While the AMBER molecular dynamics package has successfully been used for alchemical free energy simulations in academic research groups for decades, widespread impact in industrial drug discovery settings has been minimal because of the previous limitations within the AMBER alchemical code, coupled with challenges in system setup and postprocessing workflows. Through a close academia-industry collaboration we have addressed many of the previous limitations with an aim to improve accuracy, efficiency, and robustness of alchemical binding free energy simulations in industrial drug discovery applications. Here, we highlight some of the recent advances in AMBER20 with a focus on alchemical binding free energy (BFE) calculations, which are less computationally intensive than alternative binding free energy methods where full binding/unbinding paths are explored. In addition to scientific and technical advances in AMBER20, we also describe the essential practical aspects associated with running relative alchemical BFE calculations, along with recommendations for best practices, highlighting the importance not only of the alchemical simulation code but also the auxiliary functionalities and expertise required to obtain accurate and reliable results. This work is intended to provide a contemporary overview of the scientific, technical, and practical issues associated with running relative BFE simulations in AMBER20, with a focus on real-world drug discovery applications.

Solvation Thermodynamics in Different Solvents: Water-Chloroform Partition Coefficients from Grid Inhomogeneous Solvation Theory

Reliable information on partition coefficients plays a key role in drug development, as solubility decisively affects bioavailability. In a physicochemical context, the partition coefficient of a solute between two different solvents can be described as a function of solvation free energies. Hence, substantial scientific efforts have been made toward accurate predictions of solvation free energies in various solvents. The grid inhomogeneous solvation theory (GIST) facilitates the calculation of solvation free energies. In this study, we introduce an extended version of the GIST algorithm, which enables the calculation for chloroform in addition to water. Furthermore, GIST allows localization of enthalpic and entropic contributions. We test our approach by calculating partition coefficients between water and chloroform for a set of eight small molecules. We report a Pearson correlation coefficient of 0.96 between experimentally determined and calculated partition coefficients. The capability to reliably predict partition coefficients between water and chloroform and the possibility to localize their contributions allow the optimization of a compound's partition coefficient. Therefore, we presume that this methodology will be of great benefit for the efficient development of pharmaceuticals.

Biosolvation Nature of Ionic Liquids: Molecular Dynamics Simulation of Methylated Nucleobases in Hydrated 1-Ethyl-3-methylimidazolium Acetate

Solvation free energies of methylated nucleobases were calculated in pure and hydrated 1-ethyl-3-methylimidazolium acetate, [Emim][Ac], ionic liquid, and pure water using classical molecular dynamics simulations using multistate Bennett's acceptance ratio method. The calculated solvation free energies in pure water were compared with the previous experimental and theoretical findings and found to be in agreement. We observe that the solvation free energy of methylated nucleobases is more in the pure ionic liquid compared to that in the pure water and on changing the mole fraction of water in the ionic liquid, the solvation free energy decreases gradually. Comparing the Coulombic and van der Waals contribution to the solvation free energy, electrostatic contribution is more compared to that of the latter for all nucleobases. To obtain the atomistic details and explain the solvation mechanism, we calculated radial distribution functions (RDFs), spatial distribution functions (SDFs), and stacking angle distribution of cations to the nucleobases. From RDFs and SDFs, we find that the acetate anions of the ionic liquid are forming strong hydrogen bonds with the amine hydrogen atoms of the nucleobases. These hydrogen bonds contribute to the major part of the Coulombic contribution to the solvation free energy. Stacking of cations to the nucleobases is primarily due to the van der Waals contribution to the solvation free energy.

Partition coefficients of methylated DNA bases obtained from free energy calculations with molecular electron density derived atomic charges

Partition coefficients serve in various areas as pharmacology and environmental sciences to predict the hydrophobicity of different substances. Recently, they have also been used to address the accuracy of force fields for various organic compounds and specifically the methylated DNA bases. In this study, atomic charges were derived by different partitioning methods (Hirshfeld and Minimal Basis Iterative Stockholder) directly from the electron density obtained by electronic structure calculations in a vacuum, with an implicit solvation model or with explicit solvation taking the dynamics of the solute and the solvent into account. To test the ability of these charges to describe electrostatic interactions in force fields for condensed phases, the original atomic charges of the AMBER99 force field were replaced with the new atomic charges and combined with different solvent models to obtain the hydration and chloroform solvation free energies by molecular dynamics simulations. Chloroform-water partition coefficients derived from the obtained free energies were compared to experimental and previously reported values obtained with the GAFF or the AMBER-99 force field. The results show that good agreement with experimental data is obtained when the polarization of the electron density by the solvent has been taken into account, and when the energy needed to polarize the electron density of the solute has been considered in the transfer free energy. These results were further confirmed by hydration free energies of polar and aromatic amino acid side chain analogs. Comparison of the two partitioning methods, Hirshfeld-I and Minimal Basis Iterative Stockholder (MBIS), revealed some deficiencies in the Hirshfeld-I method related to the unstable isolated anionic nitrogen pro-atom used in the method. Hydration free energies and partitioning coefficients obtained with atomic charges from the MBIS partitioning method accounting for polarization by the implicit solvation model are in good agreement with the experimental values. © 2018 Wiley Periodicals, Inc.

Polarizable Multipole-Based Force Field for Aromatic Molecules and Nucleobases

Aromatic molecules with π electrons are commonly involved in chemical and biological recognitions. For example, nucleobases play central roles in DNA/RNA structure and their interactions with proteins. The delocalization of the π electrons is responsible for the high polarizability of aromatic molecules. In this work, the AMOEBA force field has been developed and applied to 5 regular nucleobases and 12 aromatic molecules. The permanent electrostatic energy is expressed as atomic multipole interactions between atom pairs, and many-body polarization is accounted for by mutually induced atomic dipoles. We have systematically investigated aromatic ring stacking and aromatic-water interactions for nucleobases and aromatic molecules, as well as base–base hydrogen-bonding pair interactions, all at various distances and orientations. van der Waals parameters were determined by comparison to the quantum mechanical interaction energy of these dimers and fine-tuned using condensed phase simulation. By comparing to quantum mechanical calculations, we show that the resulting classical potential is able to accurately describe molecular polarizability, molecular vibrational frequency, and dimer interaction energy of these aromatic systems. Condensed phase properties, including hydration free energy, liquid density, and heat of vaporization, are also in good overall agreement with experimental values. The structures of benzene liquid phase and benzene-water solution were also investigated by simulation and compared with experimental and PDB structure derived statistical results.

Molecular mechanics methods for predicting protein–ligand binding

Ligand binding affinity prediction is one of the most important applications of computational chemistry. However, accurately ranking compounds with respect to their estimated binding affinities to a biomolecular target remains highly challenging. We provide an overview of recent work using molecular mechanics energy functions to address this challenge. We briefly review methods that use molecular dynamics and Monte Carlo simulations to predict absolute and relative ligand binding free energies, as well as our own work in which we have developed a physics-based scoring method that can be applied to hundreds of thousands of compounds by invoking a number of simplifying approximations. In our previous studies, we have demonstrated that our scoring method is a promising approach for improving the discrimination between ligands that are known to bind and those that are presumed not to, in virtual screening of large compound databases. In new results presented here, we explore several improvements to our computational method including modifying the dielectric constant used for the protein and ligand interiors, and empirically scaling energy terms to compensate for deficiencies in the energy model. Future directions for further improving our physics-based scoring method are also discussed.

Theoretical model of the aqua-copper [Cu(H2O)5]+cation interactions with guanine

Pentaaqua complexes of Cu(I) with guanine were optimized at the DFT B3PW91/6-31G(d) level. For the most stable structures, vibration frequencies and NBO charges were computed followed by energy analyses. The order of individual conformers was very sensitive to the method and basis sets used for the calculation. Several conformers are practically degenerated in energy. The inclusion of an entropy term changes the order of the conformers' stability. Water molecules associated at the N9 position of guanine are favored by the inclusion of the entropy correction. Bonding energies of Cu-O(aqua) interactions were estimated to be about 60 kcal mol(-1) and for Cu-N7 bonding in the range of 75-83 kcal mol(-1). The broad range in Cu-N interaction energies demonstrates the role of induction effects caused by water molecules associated at the various sites of guanine. The charge distribution of the guanine molecule is changed remarkably by the coordination of a Cu(I) cation, which can also change the base-pairing pattern of the guanine.

Free-Energy Calculations of Protein−Ligand Cation−π and Amino−π Interactions: From Vacuum to Proteinlike Environments

To probe the role of cation-pi and amino-pi interactions in the context of protein-ligand interactions, the stability of 55 X-ray cation/amino-pi motifs involving the Ade moieties of cofactor molecules and Arg, Lys, Asn, or Gln side chains of their host protein was evaluated using quantum chemistry calculations. The conjunction of vacuum interaction energies, vibrational entropy, and solvation contributions led to identify Arg-Ade as the most favorable cation/amino-pi complex in the solvents considered, followed by Asn/Gln-Ade and Lys-Ade: their minimum interaction free energies are approximately equal to -7, -4, and -2 kcal/mol, respectively, in the solvents of dielectric constant similar to that estimated for proteins (i.e., acetone, THF, and CCl(4)). Remarkably, these free-energy values of cation/amino-pi interactions correlate well with their frequency of occurrences in protein-ligand structures, which corroborates our approach in the absence of experimental data.

Modeling DNA hydration: comparison of calculated and experimental hydration properties of nuclic acid bases

Hydration properties of individual nucleic acid bases were calculated and compared with the available experimental data. Three sets of classical potential functions (PF) used in simulations of nucleic acid hydration were juxtaposed: (i) the PF developed by Poltev and Malenkov (PM), (ii) the PF of Weiner and Kollman (WK), which together with Jorgensen's TIP3P water model are widely used in the AMBER program, and (iii) OPLS (optimized potentials for liquid simulations) developed by Jorgensen (J). The global minima of interaction energy of single water molecules with all the natural nucleic acid bases correspond to the formation of two water-base hydrogen bonds (water bridging of two hydrophilic atoms of the base). The energy values of these minima calculated via PM potentials are in somewhat better conformity with mass-spectrometric data than the values calculated via WK PF. OPLS gave much weaker water-base interactions for all compounds considered, thus these PF were not used in further computations. Monte Carlo simulations of the hydration of 9-methyladenine, 1-methyluracil and 1-methylthymine were performed in systems with 400 water molecules and periodic boundary conditions. Results of simulations with PM potentials give better agreement with experimental data on hydration energies than WK PF. Computations with PM PF of the hydration energy of keto and enol tautomers of 9-methylguanine can account for the shift in the tautomeric equilibrium of guanine in aqueous media to a dominance of the keto form in spite of nearly equal intrinsic stability of keto and enol tautomers. The results of guanine hydration computations are discussed in relation to mechanisms of base mispairing errors in nucleic acid biosynthesis. The data presented in this paper along with previous results on simulation of hydration shell structures in DNA duplex grooves provide ample evidence for the advantages of PM PF in studies of nucleic-acid hydration.

The hydration and solvent polarization effects of nucleotide bases

A combined Monte Carlo quantum mechanical and molecular mechanical (QM/MM) simulation method is used to determine the free energy of hydration and the solvent polarization effect for the nucleotide bases. In the present AM1/TIP3P model, the solute molecule is characterized by valence electrons and effective nucleus cores with Hartree-Fock molecular orbital theory incorporating a solute-solvent interaction Hamiltonian. It is found that polarization energy contributes up to 37%-61% of the total solute-solvent interaction for the systems considered. The computed free energies of hydration are compared with previous theoretical results.

Synthesis, molecular modelling, and NMR structure determination of four cyclic peptide antagonists of endothelin

A combined distance geometry and molecular mechanics/dynamics (MM/MD) protocol was unable to predict the active conformation of the cyclic pentapeptide inhibitor of endothelin-1 receptor, BQ-123, and two analogues. However, the MM/MD method alone is sufficient to predict the solution conformation of a third analogue. In that one case, the combination of proline at residue 3 and an N alpha-methyl substitution at residue 5 provides enough internal constraints to eliminate conformational flexibility seen in the other three analogues. For this constrained analogue, the 50 lowest energy conformations (out of a set of 500 DGEOM-generated, MM/MD minimized conformations) differ by no more than 3.9 kcal/mol. Thirty three of these 50 conformations have backbone atom RMSDs of less than 0.33 A, relative to the lowest energy conformation. The accuracy of this MM/MD model is verified by determining the solution structure of each of the four analogues with 2D NMR techniques. Each of the cyclic pentapeptides has a well defined solution conformation where a proline residue is clearly in a gamma-turn, leaving the remaining residues in a loose beta-turn. All four experimental NMR conformations agree closely with the MM/MD model.

Absolute and relative binding free energy calculations of the interaction of biotin and its analogs with streptavidin using molecular dynamics/free energy perturbation approaches

We present calculations of the absolute and relative binding free energies of complexation of streptavidin with biotin and its analogs by means of a thermodynamic free energy perturbation method implemented with molecular dynamics. Using the recently solved crystal structure of the streptavidin-biotin complex, biotin was mutated into a dummy molecule as well as thiobiotin and iminobiotin both in the protein and in solution. The calculated absolute binding free energy was dependent on the simulation model used. Encouragingly, the "best models" provided a reasonable semiquantitative reproduction (-20 to -22 kcal/mol) of the experimental free energy (-18.3 kcal/mol). Furthermore, the calculated results give clear insights into the binding nature of the protein-ligand complex, showing that the van der Waals energy dominates the electrostatic and hydrogen bonding energies in the binding of biotin by streptavidin. Specifically, the mutation of biotin into a dummy molecule in solution has a delta G (van der Waals) approximately -4 kcal/mol, due to the cancellation of dispersion and repulsion "cavity" effects. On the other hand, in the protein, a very small free energy price must be paid to create a cavity since one already exists and the mutation of biotin into a dummy molecule has a delta G (van der Waals) approximately 15 kcal/mol. These results are also consistent with the interpretation that the entropy increase to be expected from hydrophobic interactions from desolvation of biotin is counterbalanced by a decrease in entropy accompanying the formation of buried hydrogen bonds, which have been derived from the apparently conflicting experimental data. They provide an alternative interpretation of the reason for the extremely high affinity of the biotin-streptavidin interaction than that recently proposed by Weber et al. (J. Am. Chem. Soc. 114:3197, 1992). In the case of the relative binding free energies, the calculated values of 3.8 +/- 0.6 and 7.2 +/- 0.6 kcal/mol compare well with the experimental values of 3.6 and 6.2 kcal/mol for the perturbation of biotin to thiobiotin and iminobiotin, respectively in the related protein avidin. The calculations indicate that desolvation of the ligand is important in understanding the relative affinity of the ligands with the protein. The above successful simulations suggest that the molecular dynamics/free energy perturbation method is useful for understanding the energetic features affecting the binding between proteins and ligands, since it is generally difficult to determine these factors unambiguously by experiment.(ABSTRACT TRUNCATED AT 400 WORDS)