High-throughput Binding Affinity Calculations at Extreme Scales

Resistance to chemotherapy and molecularly targeted therapies is a major factor in limiting the effectiveness of cancer treatment. In many cases, resistance can be linked to genetic changes in target proteins, either pre-existing or evolutionarily selected during treatment. Key to overcoming this challenge is an understanding of the molecular determinants of drug binding. Using multi-stage pipelines of molecular simulations we can gain insights into the binding free energy and the residence time of a ligand, which can inform both stratified and personal treatment regimes and drug development. To support the scalable, adaptive and automated calculation of the binding free energy on high-performance computing resources, we introduce the High- throughput Binding Affinity Calculator (HTBAC). HTBAC uses a building block approach in order to attain both workflow flexibility and performance. We demonstrate close to perfect weak scaling to hundreds of concurrent multi-stage binding affinity calculation pipelines. This permits a rapid time-to-solution that is essentially invariant of the calculation protocol, size of candidate ligands and number of ensemble simulations. As such, HTBAC advances the state of the art of binding affinity calculations and protocols

Hit-to-lead and lead optimization binding free energy calculations for G protein-coupled receptors

We apply the hit-to-lead ESMACS (enhanced sampling of molecular dynamics with approximation of continuum solvent) and lead-optimization TIES (thermodynamic integration with enhanced sampling) methods to compute the binding free energies of a series of ligands at the A1 and A2A adenosine receptors, members of a subclass of the GPCR (G protein-coupled receptor) superfamily. Our predicted binding free energies, calculated using ESMACS, show a good correlation with previously reported experimental values of the ligands studied. Relative binding free energies, calculated using TIES, accurately predict experimentally determined values within a mean absolute error of approximately 1 kcal mol−1. Our methodology may be applied widely within the GPCR superfamily and to other small molecule–receptor protein systems

The performance of ensemble-based free energy protocols in computing binding affinities to ROS1 kinase

Optimization of binding affinities for compounds to their target protein is a primary objective in drug discovery. Herein we report on a collaborative study that evaluates a set of compounds binding to ROS1 kinase. We use ESMACS (enhanced sampling of molecular dynamics with approximation of continuum solvent) and TIES (thermodynamic integration with enhanced sampling) protocols to rank the binding free energies. The predicted binding free energies from ESMACS simulations show good correlations with experimental data for subsets of the compounds. Consistent binding free energy differences are generated for TIES and ESMACS. Although an unexplained overestimation exists, we obtain excellent statistical rankings across the set of compounds from the TIES protocol, with a Pearson correlation coefficient of 0.90 between calculated and experimental activities

Ensemble Simulations and Experimental Free Energy Distributions: Evaluation and Characterization of Isoxazole Amides as SMYD3 Inhibitors

Optimization of binding affinities for ligands to their target protein is a primary objective in rational drug discovery. Herein, we report on a collaborative study that evaluates various compounds designed to bind to the SET and MYND domain-containing protein 3 (SMYD3). SMYD3 is a histone methyltransferase and plays an important role in transcriptional regulation in cell proliferation, cell cycle, and human carcinogenesis. Experimental measurements using the scintillation proximity assay show that the distributions of binding free energies from a large number of independent measurements exhibit non-normal properties. We use ESMACS (enhanced sampling of molecular dynamics with approximation of continuum solvent) and TIES (thermodynamic integration with enhanced sampling) protocols to predict the binding free energies and to provide a detailed chemical insight into the nature of ligand-protein binding. Our results show that the 1-trajectory ESMACS protocol works well for the set of ligands studied here. Although one unexplained outlier exists, we obtain excellent statistical ranking across the set of compounds from the ESMACS protocol and good agreement between calculations and experiments for the relative binding free energies from the TIES protocol. ESMACS and TIES are again found to be powerful protocols for the accurate comparison of the binding free energies

Application of ESMACS binding free energy protocols to diverse datasets: Bromodomain-containing protein 4

As the application of computational methods in drug discovery pipelines becomes more widespread it is increasingly important to understand how reproducible their results are and how sensitive they are to choices made in simulation setup and analysis. Here we use ensemble simulation protocols, termed ESMACS (enhanced sampling of molecular dynamics with approximation of continuum solvent), to investigate the sensitivity of the popular molecular mechanics Poisson-Boltzmann surface area (MMPBSA) methodology. Using the bromodomain-containing protein 4 (BRD4) system bound to a diverse set of ligands as our target, we show that robust rankings can be produced only through combining ensemble sampling with multiple trajectories and enhanced solvation via an explicit ligand hydration shell

The effect of protein mutations on drug binding suggests ensuing personalised drug selection

The advent of personalised medicine promises a deeper understanding of mechanisms and therefore therapies. However, the connection between genomic sequences and clinical treatments is often unclear. We studied 50 breast cancer patients belonging to a population-cohort in the state of Qatar. From Sanger sequencing, we identified several new deleterious mutations in the estrogen receptor 1 gene (ESR1). The effect of these mutations on drug treatment in the protein target encoded by ESR1, namely the estrogen receptor, was achieved via rapid and accurate protein-ligand binding affinity interaction studies which were performed for the selected drugs and the natural ligand estrogen. Four nonsynonymous mutations in the ligand-binding domain were subjected to molecular dynamics simulation using absolute and relative binding free energy methods, leading to the ranking of the efficacy of six selected drugs for patients with the mutations. Our study shows that a personalised clinical decision system can be created by integrating an individual patient's genomic data at the molecular level within a computational pipeline which ranks the efficacy of binding of particular drugs to variant proteins

Uncertainty quantification in classical molecular dynamics

Molecular dynamics simulation is now a widespread approach for understanding complex systems on the atomistic scale. It finds applications from physics and chemistry to engineering, life and medical science. In the last decade, the approach has begun to advance from being a computer-based means of rationalizing experimental observations to producing apparently credible predictions for a number of real-world applications within industrial sectors such as advanced materials and drug discovery. However, key aspects concerning the reproducibility of the method have not kept pace with the speed of its uptake in the scientific community. Here, we present a discussion of uncertainty quantification for molecular dynamics simulation designed to endow the method with better error estimates that will enable it to be used to report actionable results. The approach adopted is a standard one in the field of uncertainty quantification, namely using ensemble methods, in which a sufficiently large number of replicas are run concurrently, from which reliable statistics can be extracted. Indeed, because molecular dynamics is intrinsically chaotic, the need to use ensemble methods is fundamental and holds regardless of the duration of the simulations performed. We discuss the approach and illustrate it in a range of applications from materials science to ligand-protein binding free energy estimation. This article is part of the theme issue 'Reliability and reproducibility in computational science: implementing verification, validation and uncertainty quantification in silico'

Enabling trade-offs between accuracy and computational cost: Adaptive algorithms to reduce time to clinical insight

The efficacy of drug treatments depends on how tightly small molecules bind to their target proteins. Quantifying the strength of these interactions (the so called 'binding affinity') is a grand challenge of computational chemistry, surmounting which could revolutionize drug design and provide the platform for patient specific medicine. Recently, evidence from blind challenge predictions and retrospective validation studies has suggested that molecular dynamics (MD) can now achieve useful predictive accuracy (1 kcal/mol) This accuracy is sufficient to greatly accelerate hit to lead and lead optimization. To translate these advances in predictive accuracy so as to impact clinical and/or industrial decision making requires that binding free energy results must be turned around on reduced timescales without loss of accuracy. This demands advances in algorithms, scalable software systems, and intelligent and efficient utilization of supercomputing resources. This work is motivated by the real world problem of providing insight from drug candidate data on a time scale that is as short as possible. Specifically, we reproduce results from a collaborative project between UCL and GlaxoSmithKline to study a congeneric series of drug candidates binding to the BRD4 protein-inhibitors of which have shown promising preclinical efficacy in pathologies ranging from cancer to inflammation. We demonstrate the use of a framework called HTBAC, designed to support the aforementioned requirements of accurate and rapid drug binding affinity calculations. HTBAC facilitates the execution of the numbers of simulations while supporting the adaptive execution of algorithms. Furthermore, HTBAC enables the selection of simulation parameters during runtime which can, in principle, optimize the use of computational resources whilst producing results within a target uncertainty

Concurrent and Adaptive Extreme Scale Binding Free Energy Calculations

The efficacy of drug treatments depends on how tightly small molecules bind to their target proteins. The rapid and accurate quantification of the strength of these interactions (as measured by binding affinity) is a grand challenge of computational chemistry, surmounting which could revolutionize drug design and provide the platform for patient-specific medicine. Recent evidence suggests that molecular dynamics (MD) can achieve useful predictive accuracy (< 1 kcal/mol). For this predictive accuracy to impact clinical decision making, binding free energy computational campaigns must provide results rapidly and without loss of accuracy. This demands advances in algorithms, scalable software systems, and efficient utilization of supercomputing resources. We introduce a framework called HTBAC, designed to support accurate and scalable drug binding affinity calculations, while marshaling large simulation campaigns. We show that HTBAC supports the specification and execution of free-energy protocols at scale. This paper makes three main contributions: (1) shows the importance of adaptive execution for ensemble-based free energy protocols to improve binding affinity accuracy; (2) presents and characterizes HTBAC -- a software system that enables the scalable and adaptive execution of binding affinity protocols at scale; and (3) for a widely used free-energy protocol (TIES), shows improvements in the accuracy of simulations for a fixed amount of resource, or reduced resource consumption for a fixed accuracy as a consequence of adaptive execution

Molecular dynamics simulation of drug resistance in HIV-1 protease and reverse transcriptase

The emergence of drug resistant strains of HIV represents a major challenge in the treatment of patients who contract the virus. We investigate the use of classical molecular dynamics to give quantitative and qualitative molecular insight into the causes of resistance in the two main drug targets in HIV, protease and reverse transcriptase. We initially establish a simulation and free energy analysis protocol for the study of resistance in protease. Focusing on the binding of the inhibitor lopinavir to a series of six mutants with increasing resistance we demonstrate that ensemble simulations exhibit significantly enhanced thermodynamic sampling over single long simulations. We achieve accurate and converged relative binding free energies, reproducible to within 0.5 kcal mol^-1. The experimentally derived ranking of the systems is reproduced with a correlation coefficient of 0.89 and a mean relative deviation from experiment of 0.9 kcal mol^-1. Our protocol is then applied to investigate a patient derived viral sequence for which contradictory resistance assessments for lopinavir were obtained from existing clinical decision support systems (CDSS). Mutations at only three locations (L10I, A71I/V and L90M) in uenced the ranking. Free energies were computed for HXB2 wildtype sequences incorporating each mutation individually and all possible combinations, along with the full patient sequence. Only in the case of the patient sequence was any resistance observed. This observation suggests an explanation for the discordance found using the CDSS. The effects on drug binding of the mutations at positions 10, 71 and 90 appear to be highly dependent on the background mutations present in the remainder of the sequence. In preparation for the extension of our simulation and free energy protocol to reverse transcriptase the impact of binding both natural DNA substrates and two non nucleoside reverse transcriptase inhibitor (NNRTI) class drugs on the dynamics of reverse transcriptase are investigated. Free energies of both inhibitors (efavirenz and neviripine) are determined which are seen to be independent of the subdomain motions of the protein observed during simulation. Preliminary calculations of the free energies for a set of NNRTI resistant mutants bound to efavirenz are also presented