Computing Clinically Relevant Binding Free Energies of HIV-1 Protease Inhibitors

Abstract

The use of molecular simulation to estimate the strength of macromolecular binding free energies is becoming increasingly widespread, with goals ranging from lead optimization and enrichment in drug discovery to personalizing or stratifying treatment regimes. In order to realize the potential of such approaches to predict new results, not merely to explain previous experimental findings, it is necessary that the methods used are reliable and accurate, and that their limitations are thoroughly understood. However, the computational cost of atomistic simulation techniques such as molecular dynamics (MD) has meant that until recently little work has focused on validating and verifying the available free energy methodologies, with the consequence that many of the results published in the literature are not reproducible. Here, we present a detailed analysis of two of the most popular approximate methods for calculating binding free energies from molecular simulations, molecular mechanics Poisson–Boltzmann surface area (MMPBSA) and molecular mechanics generalized Born surface area (MMGBSA), applied to the nine FDA-approved HIV-1 protease inhibitors. Our results show that the values obtained from replica simulations of the same protease–drug complex, differing only in initially assigned atom velocities, can vary by as much as 10 kcal mol–1, which is greater than the difference between the best and worst binding inhibitors under investigation. Despite this, analysis of ensembles of simulations producing 50 trajectories of 4 ns duration leads to well converged free energy estimates. For seven inhibitors, we find that with correctly converged normal mode estimates of the configurational entropy, we can correctly distinguish inhibitors in agreement with experimental data for both the MMPBSA and MMGBSA methods and thus have the ability to rank the efficacy of binding of this selection of drugs to the protease (no account is made for free energy penalties associated with protein distortion leading to the over estimation of the binding strength of the two largest inhibitors ritonavir and atazanavir). We obtain improved rankings and estimates of the relative binding strengths of the drugs by using a novel combination of MMPBSA/MMGBSA with normal mode entropy estimates and the free energy of association calculated directly from simulation trajectories. Our work provides a thorough assessment of what is required to produce converged and hence reliable free energies for protein–ligand binding.

Novel chemical scaffolds of the tumor marker AKR1B10 inhibitors discovered by 3D QSAR pharmacophore modeling

“…The MM-PBSA analyses of the binding energies showed lower energies for the AKR1B10-ligand complexes using hit 1, hit 2, and hit 3 as the ligand and the highest binding affinities due to the presence of hydrogen bonding and hydrophobic interactions in each of these complexes. The free energy calculations using MM/PBSA method have shown to produce reliable and reproducible predictions of binding affinity that could be used in drug discovery and personalized medicine applications [48] . It can be helpful to rank the lead compounds according to their predicted affinities to the target and this information may be used to guide molecular design and synthesis [49] .…”

Rapid, Precise, and Reproducible Prediction of Peptide–MHC Binding Affinities from Molecular Dynamics That Correlate Well with Experiment

“…A sufficient number of replicas is needed to ensue full error control and hence reproducibility. 50 replica ensembles are used in this study, the same size we have used in our previous studies, although the study on HIV-1 protease inhibitors indicated that convergence can be reached for ensembles with a smaller number of replicas 14 One particularly interesting feature of the ESMACS protocol we report on here is the enhancement of predictions when we move from the so-called 1-trajectory to the 3-trajectory variant of the method, relaxing the assumption that the substrate and receptor protein sample similar conformations in both the free and bound states. The 1-trajectory method, in which the energies of protein-inhibitor complexes, receptor proteins and ligands/drugs/peptides are extracted from the MD trajectories of the complexes alone, has been used in virtually all MMPBSA studies in the literature.…”

“…Our approach is based on a 'protocol' originally reported for estimating binding affinities of inhibitors to HIV-1 protease 14,19 . We call this protocol "enhanced sampling of molecular 7 dynamics with the approximation of continuum solvent" (ESMACS).…”

Stability of Iowa mutant and wild type Aβ-peptide aggregates

“…46,47 In order to gain further insight into the relative stability of our four oligomers, we have done MM-PBSA calculations of the oligomers allowing us to monitor their interactions 48 through calculating approximate binding free energies from molecular simulations. 49,50 While the MM-PBSA approach in general does not replicate the absolute binding free energy values, 51 we chose this approach because it allows one to calculate quickly an estimate for differences in the free energy of binding, and because it usually exhibits a good correlation with experimental data. 52 In the present study, the binding energy between the two β-sheets (that is between the pentamers that form the decamer) is estimated with the MM-PBSA methodology as implemented in AMBER12.…”