Accurate and Reliable Prediction of the Binding Affinities of Macrocycles to Their Protein Targets

Identification of Novel Potent NSD2-PWWP1 Ligands Using Structure-Based Design and Computational Approaches

Dysregulation of histone methyl transferase nuclear receptor-binding SET domain 2 (NSD2) has been implicated in several hematological and solid malignancies. NSD2 is a large multidomain protein that carries histone writing and histone reading functions. To date, identifying inhibitors of the enzymatic activity of NSD2 has proven challenging in terms of potency and SET domain selectivity. Inhibition of the NSD2-PWWP1 domain using small molecules has been considered as an alternative approach to reduce NSD2-unregulated activity. In this article, we present novel computational chemistry approaches, encompassing free energy perturbation coupled to machine learning (FEP/ML) models as well as virtual screening (VS) activities, to identify high-affinity NSD2 PWWP1 binders. Through these activities, we have identified the most potent NSD2-PWWP1 binder reported so far in the literature: compound 34 (pIC50 = 8.2). The compounds identified herein represent useful tools for studying the role of PWWP1 domains for inhibition of human NSD2.

Selection pressures on evolution of ribonuclease H explored with rigorous free-energy-based design

Understanding natural protein evolution and designing novel proteins are motivating interest in development of high-throughput methods to explore large sequence spaces. In this work, we demonstrate the application of multisite λ dynamics (MSλD), a rigorous free energy simulation method, and chemical denaturation experiments to quantify evolutionary selection pressure from sequence-stability relationships and to address questions of design. This study examines a mesophilic phylogenetic clade of ribonuclease H (RNase H), furthering its extensive characterization in earlier studies, focusing on E. coli RNase H (ecRNH) and a more stable consensus sequence (AncCcons) differing at 15 positions. The stabilities of 32,768 chimeras between these two sequences were computed using the MSλD framework. The most stable and least stable chimeras were predicted and tested along with several other sequences, revealing a designed chimera with approximately the same stability increase as AncCcons, but requiring only half the mutations. Comparing the computed stabilities with experiment for 12 sequences reveals a Pearson correlation of 0.86 and root mean squared error of 1.18 kcal/mol, an unprecedented level of accuracy well beyond less rigorous computational design methods. We then quantified selection pressure using a simple evolutionary model in which sequences are selected according to the Boltzmann factor of their stability. Selection temperatures from 110 to 168 K are estimated in three ways by comparing experimental and computational results to evolutionary models. These estimates indicate selection pressure is high, which has implications for evolutionary dynamics and for the accuracy required for design, and suggests accurate high-throughput computational methods like MSλD may enable more effective protein design.

The maximal and current accuracy of rigorous protein-ligand binding free energy calculations

Computational techniques can speed up the identification of hits and accelerate the development of candidate molecules for drug discovery. Among techniques for predicting relative binding affinities, the most consistently accurate is free energy perturbation (FEP), a class of rigorous physics-based methods. However, uncertainty remains about how accurate FEP is and can ever be. Here, we present what we believe to be the largest publicly available dataset of proteins and congeneric series of small molecules, and assess the accuracy of the leading FEP workflow. To ascertain the limit of achievable accuracy, we also survey the reproducibility of experimental relative affinity measurements. We find a wide variability in experimental accuracy and a correspondence between binding and functional assays. When careful preparation of protein and ligand structures is undertaken, FEP can achieve accuracy comparable to experimental reproducibility. Throughout, we highlight reliable protocols that can help maximize the accuracy of FEP in prospective studies.

An Intriguing Purview on the Design of Macrocyclic Inhibitors for Unexplored Protein Kinases through Their Binding Site Comparison

Kinases play an important role in regulating various intracellular signaling pathways that control cell proliferation, differentiation, survival, and other cellular processes, and their deregulation causes more than 400 diseases. Consequently, macrocyclization can be considered a noteworthy approach to developing new therapeutic agents for human diseases. Macrocyclization has emerged as an effective drug discovery strategy over the past decade to improve target selectivity and potency of small molecules. Small compounds with linear structures upon macrocyclization can lead to changes in their physicochemical and biological properties by firmly reducing conformational flexibility. A number of distinct protein kinases exhibit similar binding sites. Comparison of protein binding sites provides crucial insights for drug discovery and development. Binding site similarities are helpful in understanding polypharmacology, identifying potential off-targets, and repurposing known drugs. In this review, we focused on comparing the binding sites of those kinases for which macrocyclic inhibitors are available/studied so far. Furthermore, we calculated the volume of the binding site pocket for each targeted kinase and then compared it with the binding site pocket of the kinase for which only acyclic inhibitors were designed to date. Our review and analysis of several explored kinases might be useful in targeting new protein kinases for macrocyclic drug discovery.

Mycobacterium tuberculosis PptT Inhibitors Based on Heterocyclic Replacements of Amidinoureas

4'-Phosphopantetheinyl transferase (PptT) is an essential enzyme for Mycobacterium tuberculosis (Mtb) survival and virulence and therefore an attractive target for a tuberculosis therapeutic. In this work, two modeling-informed approaches toward the isosteric replacement of the amidinourea moiety present in the previously reported PptT inhibitor AU 8918 are reported. Although a designed 3,5-diamino imidazole unexpectedly adopted an undesired tautomeric form and was inactive, replacement of the amidinourea moiety afforded a series of active PptT inhibitors containing 2,6-diaminopyridine scaffolds.

Enhancing Hit Discovery in Virtual Screening through Absolute Protein-Ligand Binding Free-Energy Calculations

In the hit identification stage of drug discovery, a diverse chemical space needs to be explored to identify initial hits. Contrary to empirical scoring functions, absolute protein-ligand binding free-energy perturbation (ABFEP) provides a theoretically more rigorous and accurate description of protein-ligand binding thermodynamics and could, in principle, greatly improve the hit rates in virtual screening. In this work, we describe an implementation of an accurate and reliable ABFEP method in FEP+. We validated the ABFEP method on eight congeneric compound series binding to eight protein receptors including both neutral and charged ligands. For ligands with net charges, the alchemical ion approach is adopted to avoid artifacts in electrostatic potential energy calculations. The calculated binding free energies correlate with experimental results with a weighted average of R² = 0.55 for the entire dataset. We also observe an overall root-mean-square error (RMSE) of 1.1 kcal/mol after shifting the zero-point of the simulation data to match the average experimental values. Through ABFEP calculations using apo versus holo protein structures, we demonstrated that the protein conformational and protonation state changes between the apo and holo proteins are the main physical factors contributing to the protein reorganization free energy manifested by the overestimation of raw ABFEP calculated binding free energies using the holo structures of the proteins. Furthermore, we performed ABFEP calculations in three virtual screening applications for hit enrichment. ABFEP greatly improves the hit rates as compared to docking scores or other methods like metadynamics. The good performance of ABFEP in rank ordering compounds demonstrated in this work confirms it as a useful tool to improve the hit rates in virtual screening, thus facilitating hit discovery.

Prioritizing Small Sets of Molecules for Synthesis through in-silico Tools: A Comparison of Common Ranking Methods

Prioritizing molecules for synthesis is a key role of computational methods within medicinal chemistry. Multiple tools exist for ranking molecules, from the cheap and popular molecular docking methods to more computationally expensive molecular-dynamics (MD)-based methods. It is often questioned whether the accuracy of the more rigorous methods justifies the higher computational cost and associated calculation time. Here, we compared the performance on ranking the binding of small molecules for seven scoring functions from five docking programs, one end-point method (MM/GBSA), and two MD-based free energy methods (PMX, FEP+). We investigated 16 pharmaceutically relevant targets with a total of 423 known binders. The performance of docking methods for ligand ranking was strongly system dependent. We observed that MD-based methods predominantly outperformed docking algorithms and MM/GBSA calculations. Based on our results, we recommend the application of MD-based free energy methods for prioritization of molecules for synthesis in lead optimization, whenever feasible.

Best practices for constructing, preparing, and evaluating protein-ligand binding affinity benchmarks [Article v0.1]

Free energy calculations are rapidly becoming indispensable in structure-enabled drug discovery programs. As new methods, force fields, and implementations are developed, assessing their expected accuracy on real-world systems (benchmarking) becomes critical to provide users with an assessment of the accuracy expected when these methods are applied within their domain of applicability, and developers with a way to assess the expected impact of new methodologies. These assessments require construction of a benchmark-a set of well-prepared, high quality systems with corresponding experimental measurements designed to ensure the resulting calculations provide a realistic assessment of expected performance when these methods are deployed within their domains of applicability. To date, the community has not yet adopted a common standardized benchmark, and existing benchmark reports suffer from a myriad of issues, including poor data quality, limited statistical power, and statistically deficient analyses, all of which can conspire to produce benchmarks that are poorly predictive of real-world performance. Here, we address these issues by presenting guidelines for (1) curating experimental data to develop meaningful benchmark sets, (2) preparing benchmark inputs according to best practices to facilitate widespread adoption, and (3) analysis of the resulting predictions to enable statistically meaningful comparisons among methods and force fields. We highlight challenges and open questions that remain to be solved in these areas, as well as recommendations for the collection of new datasets that might optimally serve to measure progress as methods become systematically more reliable. Finally, we provide a curated, versioned, open, standardized benchmark set adherent to these standards (PLBenchmarks) and an open source toolkit for implementing standardized best practices assessments (arsenic) for the community to use as a standardized assessment tool. While our main focus is free energy methods based on molecular simulations, these guidelines should prove useful for assessment of the rapidly growing field of machine learning methods for affinity prediction as well.

RestraintMaker: a graph-based approach to select distance restraints in free-energy calculations with dual topology

The calculation of relative binding free energies (RBFE) involves the choice of the end-state/system representation, of a sampling approach, and of a free-energy estimator. System representations are usually termed “single topology” or “dual topology”. As the terminology is often used ambiguously in the literature, a systematic categorization of the system representations is proposed here. In the dual-topology approach, the molecules are simulated as separate molecules. Such an approach is relatively easy to automate for high-throughput RBFE calculations compared to the single-topology approach. Distance restraints are commonly applied to prevent the molecules from drifting apart, thereby improving the sampling efficiency. In this study, we introduce the program RestraintMaker, which relies on a greedy algorithm to find (locally) optimal distance restraints between pairs of atoms based on geometric measures. The algorithm is further extended for multi-state methods such as enveloping distribution sampling (EDS) or multi-site \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda $$\end{document}λ-dynamics. The performance of RestraintMaker is demonstrated for toy models and for the calculation of relative hydration free energies. The Python program can be used in script form or through an interactive GUI within PyMol. The selected distance restraints can be written out in GROMOS or GROMACS file formats. Additionally, the program provides a human-readable JSON format that can easily be parsed and processed further. The code of RestraintMaker is freely available on GitHub https://github.com/rinikerlab/restraintmaker.

Computational studies of molecular pre-organization through macrocyclization: Conformational distribution analysis of closely related non-macrocyclic and macrocyclic analogs

Macrocycles form an important compound class in medicinal chemistry due to their interesting structural and biological properties. To help design macrocycles, it is important to understand how the conformational preferences are affected upon macrocyclization of a lead compound. To address this, we collected a unique data set of protein-ligand complexes containing "non-macrocyclic" ("linear") ligands matched with macrocyclic analogs binding to the same protein in a similar pose. Out of the 39 co-crystallized ligands considered, 10 were linear and 29 were macrocyclic. To enable a more general analysis, 128 additional ligands from the publications associated with these protein data bank entries were added to the data set. Using in total 167 collected ligands, we investigated if the conformers in the macrocyclic conformational ensembles were more similar to the bioactive conformation in comparison to the conformers of their linear counterparts. Unexpectedly, in most cases the macrocycle conformational ensemble distributions were not very different from those of the linear compounds. Thus, care should be taken when designing macrocycles with the aim to focus their conformational preference towards the bioactive conformation. We also set out to investigate potential conformational flexibility differences between the two compound classes, computational energy window settings and evaluate a literature metric for approximating the conformational focusing on the bioactive conformation.

Optimization of New Catalytic Topoisomerase II Inhibitors as an Anti-Cancer Therapy

Simple Summary

DNA topoisomerase II (TOP2) is a drug target for many types of cancers. However, clinically used TOP2 inhibitors not only kill cancer cells, but also damage normal cells, and can even give rise to other types of cancers. To discover new TOP2 inhibitors to more effectively treat cancer patients, we have applied computer-aided drug design technology to develop several TOP2 inhibitors that can strongly inhibit cancer cell growth but exert low side effects. Results of one exemplary compound are presented in this study. It shows several promising drug-like properties that can be potentially developed into anticancer drugs.

Abstract

Clinically used topoisomerase II (TOP2) inhibitors are poison inhibitors that induce DNA damage to cause cancer cell death. However, they can also destroy benign cells and thereby show serious side effects, including cardiotoxicity and drug-induced secondary malignancy. New TOP2 inhibitors with a different mechanism of action (MOA), such as catalytic TOP2 inhibitors, are needed to more effectively control tumor growth. We have applied computer-aided drug design to develop a new group of small molecule inhibitors that are derivatives of our previously identified lead compound T60. Particularly, the compound T638 has shown improved solubility and microsomal stability. It is a catalytic TOP2 inhibitor that potently suppresses TOP2 activity. T638 has a novel MOA by which it binds TOP2 proteins and blocks TOP2–DNA interaction. T638 strongly inhibits cancer cell growth, but exhibits limited genotoxicity to cells. These results indicate that T638 is a promising drug candidate that warrants further development into clinically used anticancer drugs.

Scaffold Hopping Transformations Using Auxiliary Restraints for Calculating Accurate Relative Binding Free Energies

In silico screening of drug-target interactions is a key part of the drug discovery process. Changes in the drug scaffold via contraction or expansion of rings, the breaking of rings, and the introduction of cyclic structures from acyclic structures are commonly applied by medicinal chemists to improve binding affinity and enhance favorable properties of candidate compounds. These processes, commonly referred to as scaffold hopping, are challenging to model computationally. Although relative binding free energy (RBFE) calculations have shown success in predicting binding affinity changes caused by perturbing R-groups attached to a common scaffold, applications of RBFE calculations to modeling scaffold hopping are relatively limited. Scaffold hopping inevitably involves breaking and forming bond interactions of quadratic functional forms, which is highly challenging. A novel method for handling ring opening/closure/contraction/expansion and linker contraction/expansion is presented here. To the best of our knowledge, RBFE calculations on linker contraction/expansion have not been previously reported. The method uses auxiliary restraints to hold the atoms at the ends of a bond in place during the breaking and forming of the bonds. The broad applicability of the method was demonstrated by examining perturbations involving small-molecule macrocycles and mutations of proline in proteins. High accuracy was obtained using the method for most of the perturbations studied. The rigor of the method was isolated from the force field by validating the method using relative and absolute hydration free energy calculations compared to standard simulation results. Unlike other methods that rely on λ-dependent functional forms for bond interactions, the method presented here can be employed using modern molecular dynamics software without modification of codes or force field functions.

A strategy for proline and glycine mutations to proteins with alchemical free energy calculations

Computation of the thermodynamic consequences of protein mutations holds great promise in protein biophysics and design. Alchemical free energy methods can give improved estimates of mutational free energies, and are already widely used in calculations of relative and absolute binding free energies in small molecule design problems. In principle, alchemical methods can address any amino acid mutation with an appropriate alchemical pathway, but identifying a strategy that produces such a path for proline and glycine mutations is an ongoing challenge. Most current strategies perturb only side chain atoms, while proline and glycine mutations also alter the backbone parameters and backbone ring topology. Some strategies also perturb backbone parameters and enable glycine mutations. This work presents a strategy that enables both proline and glycine mutations and comprises two key elements: a dual backbone with restraints and scaling of bonded terms, facilitating backbone parameter changes, and a soft bond in the proline ring, enabling ring topology changes in proline mutations. These elements also have utility for core hopping and macrocycle studies in computer-aided drug design. This new strategy shows slight improvements over an alternative side chain perturbation strategy for a set T4 lysozyme mutations lacking proline and glycine, and yields good agreement with experiment for a set of T4 lysozyme proline and glycine mutations not previously studied. To our knowledge this is the first report comparing alchemical predictions of proline mutations with experiment. With this strategy in hand, alchemical methods now have access to the full palette of amino acid mutations.

CHARMM-GUI Free Energy Calculator for Absolute and Relative Ligand Solvation and Binding Free Energy Simulations

Alchemical free energy simulations have long been utilized to predict free energy changes for binding affinity and solubility of small molecules. However, while the theoretical foundation of these methods is well established, seamlessly handling many of the practical aspects regarding the preparation of the different thermodynamic end states of complex molecular systems and the numerous processing scripts often remains a burden for successful applications. In this work, we present CHARMM-GUI Free Energy Calculator (http://www.charmm-gui.org/input/fec) that provides various alchemical free energy perturbation molecular dynamics (FEP/MD) systems with input and post-processing scripts for NAMD and GENESIS. Four submodules are available: Absolute Ligand Binder (for absolute ligand binding FEP/MD), Relative Ligand Binder (for relative ligand binding FEP/MD), Absolute Ligand Solvator (for absolute ligand solvation FEP/MD), and Relative Ligand Solvator (for relative ligand solvation FEP/MD). Each module is designed to build multiple systems of a set of selected ligands at once for high-throughput FEP/MD simulations. The capability of Free Energy Calculator is illustrated by absolute and relative solvation FEP/MD of a set of ligands and absolute and relative binding FEP/MD of a set of ligands for T4-lysozyme in solution and the adenosine A_2A receptor in a membrane. The calculated free energy values are overall consistent with the experimental and published free energy results (within ∼1 kcal/mol). We hope that Free Energy Calculator is useful to carry out high-throughput FEP/MD simulations in the field of biomolecular sciences and drug discovery.

Automated Design of Macrocycles for Therapeutic Applications: From Small Molecules to Peptides and Proteins

Macrocycles and cyclic peptides are increasingly attractive therapeutic modalities as they often have improved affinity, are able to bind to extended protein surfaces, and otherwise have favorable properties. Macrocyclization of a known binder may stabilize its bioactive conformation and improve its metabolic stability, cell permeability, and in certain cases oral bioavailability. Herein, we present implementation and application of an approach that automatically generates, evaluates, and proposes cyclizations utilizing a library of well-established chemical reactions and reagents. Using the three-dimensional (3D) conformation of the linear molecule in complex with a target protein as the starting point, this approach identifies attachment points, generates linkers, evaluates their geometric compatibility, and ranks the resulting molecules with respect to their predicted conformational stability and interactions with the target protein. As we show here with prospective and retrospective case studies, this procedure can be applied for the macrocyclization of small molecules and peptides and even PROteolysis TArgeting Chimeras (PROTACs) and proteins.

Estimation of the ligand-binding free energy of checkpoint kinase 1 via non-equilibrium MD simulations

Checkpoint kinase 1 (CHK1) is a serine/threonine-protein kinase that is involved in cell cycle regulation in eukaryotes. Inhibition of CHK1 is thus considered as a promising approach in cancer therapy. In this study, the fast pulling of ligand (FPL) process was applied to predict the relative binding affinities of CHK1 inhibitors using non-equilibrium molecular dynamics (MD) simulations. The work of external harmonic forces to pull the ligand out of the binding cavity strongly correlated with the experimental binding affinity of CHK1 inhibitors with the correlation coefficient of R = -0.88 and an overall root mean square error (RMSE) of 0.99 kcal/mol. The data indicate that the FPL method is highly accurate in predicting the relative binding free energies of CHK1 inhibitors with an affordable CPU time. A new set of molecules were designed based on the molecular modeling of interactions between the known inhibitor and CHK1 as inhibitory candidates. Molecular docking and FPL results exhibited that the binding affinities of developed ligands were similar to the known inhibitor in interaction with the catalytic site of CHK1, producing very potential CHK1 inhibitors of that the inhibitory activities should be further evaluated in vitro.

Improving the Accuracy of Protein Thermostability Predictions for Single Point Mutations

Accurately predicting the protein thermostability changes upon single point mutations in silico is a challenge that has implications for understanding diseases as well as industrial applications of protein engineering. Free energy perturbation (FEP) has been applied to predict the effect of single point mutations on protein stability for over 40 years and emerged as a potentially reliable prediction method with reasonable throughput. However, applications of FEP in protein stability calculations in industrial settings have been hindered by a number of limitations, including the inability to model mutations to and from prolines in which the bonded topology of the backbone is modified and the complexity in modeling charge-changing mutations. In this study, we have extended the FEP+ protocol to enable the accurate modeling of the effects on protein stability from proline mutations and from charge-changing mutations. We also evaluated the influence of the unfolded model in the stability calculations using increasingly longer peptides with native sequence and conformations. With the abovementioned improvements, the accuracy of FEP predictions of protein stability over a data set of 87 mutations on five different proteins has drastically improved compared with previous studies, with a mean unsigned error of 0.86 kcal/mol and root mean square error of 1.11 kcal/mol, comparable with the accuracy of previously published state-of-the-art small-molecule relative binding affinity calculations, which have been shown to be capable of driving discovery projects.

Adapting free energy perturbation simulations for large macrocyclic ligands: how to dissect contributions from direct binding and free ligand flexibility

Large and flexible ligands gain increasing interest in the development of bioactive agents. They challenge the applicability of computational ligand optimization strategies originally developed for small molecules. Free energy perturbation (FEP) is often used for predicting binding affinities of small molecule ligands, however, its use for more complex ligands remains limited. Herein, we report the structure-based design of peptide macrocycles targeting the protein binding site of human adaptor protein 14-3-3. We observe a surprisingly strong dependency of binding affinities on relatively small variations in substituent size. FEP was performed to rationalize observed trends. To account for insufficient convergence of FEP, restrained calculations were performed and complemented with extensive REST MD simulations of the free ligands. These calculations revealed that changes in affinity originate both from altered direct interactions and conformational changes of the free ligand. In addition, MD simulations provided the basis to rationalize unexpected trends in ligand lipophilicity. We also verified the anticipated interaction site and binding mode for one of the high affinity ligands by X-ray crystallography. The introduced fully-atomistic simulation protocol can be used to rationalize the development of structurally complex ligands which will support future ligand maturation efforts.

D3R Grand Challenge 4: ligand similarity and MM-GBSA-based pose prediction and affinity ranking for BACE-1 inhibitors

The Drug Design Data Resource (D3R) Grand Challenges present an opportunity to assess, in the context of a blind predictive challenge, the accuracy and the limits of tools and methodologies designed to help guide pharmaceutical drug discovery projects. Here, we report the results of our participation in the D3R Grand Challenge 4 (GC4), which focused on predicting the binding poses and affinity ranking for compounds targeting the [Formula: see text]-amyloid precursor protein (BACE-1). Our ligand similarity-based protocol using HYBRID (OpenEye Scientific Software) successfully identified poses close to the native binding mode for most of the ligands with less than 2 Å RMSD accuracy. Furthermore, we compared the performance of our HYBRID-based approach to that of AutoDock Vina and DOCK 6 and found that using a reference ligand to guide the docking process is a better strategy for pose prediction and helped HYBRID to perform better here. We also conducted end-point free energy estimates on molecules dynamics based ensembles of protein-ligand complexes using molecular mechanics combined with generalized Born surface area method (MM-GBSA). We found that the binding affinity ranking based on MM-GBSA scores have poor correlation with the experimental values. Finally, the main lessons from our participation in D3R GC4 are: (i) the generation of the macrocyclic conformers is a key step for successful pose prediction, (ii) the protonation states of the BACE-1 binding site should be treated carefully, (iii) the MM-GBSA method could not discriminate well between different predicted binding poses, and (iv) the MM-GBSA method does not perform well at predicting protein-ligand binding affinities here.

D3R grand challenge 4: blind prediction of protein-ligand poses, affinity rankings, and relative binding free energies

The Drug Design Data Resource (D3R) aims to identify best practice methods for computer aided drug design through blinded ligand pose prediction and affinity challenges. Herein, we report on the results of Grand Challenge 4 (GC4). GC4 focused on proteins beta secretase 1 and Cathepsin S, and was run in an analogous manner to prior challenges. In Stage 1, participant ability to predict the pose and affinity of BACE1 ligands were assessed. Following the completion of Stage 1, all BACE1 co-crystal structures were released, and Stage 2 tested affinity rankings with co-crystal structures. We provide an analysis of the results and discuss insights into determined best practice methods.

Best Practices for Alchemical Free Energy Calculations [Article v1.0]

Alchemical free energy calculations are a useful tool for predicting free energy differences associated with the transfer of molecules from one environment to another. The hallmark of these methods is the use of "bridging" potential energy functions representing alchemical intermediate states that cannot exist as real chemical species. The data collected from these bridging alchemical thermodynamic states allows the efficient computation of transfer free energies (or differences in transfer free energies) with orders of magnitude less simulation time than simulating the transfer process directly. While these methods are highly flexible, care must be taken in avoiding common pitfalls to ensure that computed free energy differences can be robust and reproducible for the chosen force field, and that appropriate corrections are included to permit direct comparison with experimental data. In this paper, we review current best practices for several popular application domains of alchemical free energy calculations performed with equilibrium simulations, in particular relative and absolute small molecule binding free energy calculations to biomolecular targets.

Approaching protein design with multisite λ dynamics: Accurate and scalable mutational folding free energies in T4 lysozyme

The estimation of changes in free energy upon mutation is central to the problem of protein design. Modern protein design methods have had remarkable success over a wide range of design targets, but are reaching their limits in ligand binding and enzyme design due to insufficient accuracy in mutational free energies. Alchemical free energy calculations have the potential to supplement modern design methods through more accurate molecular dynamics based prediction of free energy changes, but suffer from high computational cost. Multisite λ dynamics (MSλD) is a particularly efficient and scalable free energy method with potential to explore combinatorially large sequence spaces inaccessible with other free energy methods. This work aims to quantify the accuracy of MSλD and demonstrate its scalability. We apply MSλD to the classic problem of calculating folding free energies in T4 lysozyme, a system with a wealth of experimental measurements. Single site mutants considering 32 mutations show remarkable agreement with experiment with a Pearson correlation of 0.914 and mean unsigned error of 1.19 kcal/mol. Multisite mutants in systems with up to five concurrent mutations spanning 240 different sequences show comparable agreement with experiment. These results demonstrate the promise of MSλD in exploring large sequence spaces for protein design.

Blinded prediction of protein-ligand binding affinity using Amber thermodynamic integration for the 2018 D3R grand challenge 4

In the framework of the 2018 Drug Design Data Resource grand challenge 4, blinded predictions on relative binding free energy were performed for a set of 39 ligands of the Cathepsin S protein. We leveraged the GPU-accelerated thermodynamic integration of Amber 18 to advance our computational prediction. When our entry was compared to experimental results, a good correlation was observed (Kendall's τ: 0.62, Spearman's ρ: 0.80 and Pearson's R: 0.82). We designed a parallelized transformation map that placed ligands into several groups based on common alchemical substructures; TI transformations were carried out for each ligand to the relevant substructure, and between substructures. Our calculations were all conducted using the linear potential scaling scheme in Amber TI because we believe the softcore potential/dual-topology approach as implemented in current Amber TI is highly fault-prone for some transformations. The issue is illustrated by using two examples in which typical preparation for the dual-topology approach of Amber TI fails. Overall, the high accuracy of our prediction is a result of recent advances in force fields (ff14SB and GAFF), as well as rapid calculation of ensemble averages enabled by the GPU implementation of Amber. The success shown here in a blinded prediction strongly suggests that alchemical free energy calculation in Amber is a promising tool for future commercial drug design.

Free Ligand 1D NMR Conformational Signatures To Enhance Structure Based Drug Design of a Mcl-1 Inhibitor (AZD5991) and Other Synthetic Macrocycles

The three-dimensional conformations adopted by a free ligand in solution impact bioactivity and physicochemical properties. Solution 1D NMR spectra inherently contain information on ligand conformational flexibility and three-dimensional shape, as well as the propensity of the free ligand to fully preorganize into the bioactive conformation. Herein we discuss some key learnings, distilled from our experience developing potent and selective synthetic macrocyclic inhibitors, including Mcl-1 clinical candidate AZD5991. Case studies have been selected from recent oncology research projects, demonstrating how 1D NMR conformational signatures can complement X-ray protein-ligand structural information to guide medicinal chemistry optimization. Learning to extract free ligand conformational information from routinely available 1D NMR signatures has proven to be fast enough to guide medicinal chemistry decisions within design cycles for compound optimization.

Application of Free Energy Perturbation (FEP+) to Understanding Ligand Selectivity: A Case Study to Assess Selectivity Between Pairs of Phosphodiesterases (PDE's)

Cyclic nucleotide phosphodiesterases (PDE's) are metalloenzymes that play a key role in regulating the levels of the ubiquitous second messengers, cyclic adenosine monophosphate (cAMP) and cyclic guanosine monophosphate (cGMP). In humans, 11 PDE protein families mediate numerous biochemical pathways throughout the body and are effective drug targets for the treatment of diseases ranging from central nervous system disorders to heart and pulmonary diseases. PDE's also share a highly conserved catalytic site (about 50%), thus making the design of selective drug candidates very challenging with classical structure-based design approaches given also the lack of publicly available co-crystal structures of pairs of PDE's with consistent biological data to be compared, as we show in our work. In this retrospective study, we apply free energy perturbation (FEP+) to predict the selectivity of inhibitors that bind two pairs of closely related PDE families: PDE9/1 and PDE5/6 where only 1 co-crystal structure per pair is publicly available. As another challenge, the p K_a of the PDE5/6 inhibitor is close to the experimental pH, making unclear the exact protonation state that should be used in the computational workflow. We demonstrate that running FEP+ on homology models constructed for these metalloenzymes accurately reproduces experimentally observed selectivity profiles also addressing the unclear protonation state to be used during computation with our recently developed p K_a-correction method. Based on these data, we conclude that FEP+ is a robust method for prediction of selectivity for this target class and may be helpful to address related lead optimization challenges in drug discovery.

Understanding PIM-1 kinase inhibitor interactions with free energy simulation

The proviral integration site of the Moloney leukemia virus (PIM) family includes three homologous members. PIM-1 kinase is an important target in effective therapeutic interventions of lymphomas, prostate cancer and leukemia. In the current work, we performed free energy calculations to calculate the binding affinities of several inhibitors targeting this protein. The alchemical method with integration and perturbation-based estimators and the end-point methods were compared. The computational results indicated that the alchemical method can accurately predict the binding affinities, while the end-point methods give relatively unreliable predictions. Decomposing the free energy difference into enthalpic and entropic components with MBAR reweighting enabled us to investigate the detailed thermodynamic parameters with which the entropy-enthalpy compensation in this protein-ligand binding case is identified. We then studied the conformational ensemble, and the important protein-ligand interactions were identified. The current work sheds light on the understanding of the PIM-1-kinase-inhibitor interactions at the atomic level and will be useful in the further development of potential drugs.

Protein-Ligand Binding Free Energy Calculations with FEP

Accurate and reliable calculation of protein-ligand binding free energy is of central importance in computational biophysics and structure-based drug design. Among the various methods to calculate protein-ligand binding affinities, alchemical free energy perturbation (FEP) calculations performed by way of explicitly solvated molecular dynamics simulations (FEP/MD) provide a thermodynamically rigorous and complete description of the binding event and should in turn yield highly accurate predictions. Although the original theory of FEP was proposed more than 60 years ago, subsequent applications of FEP to compute protein-ligand binding free energies in the context of drug discovery projects over much of that time period was sporadic and generally unsuccessful. This was mainly due to the limited accuracy of the available force fields, inadequate sampling of the protein-ligand conformational space, complexity of simulation set up and analysis, and the large computational resources required to pursue such calculations. Over the past few years, there have been advances in computing power, classical force field accuracy, enhanced sampling algorithms, and simulation setup. This has led to newer FEP implementations such as the FEP+ technology developed by Schrödinger Inc., which has enabled accurate and reliable calculations of protein-ligand binding free energies and positioned free energy calculations to play a guiding role in small-molecule drug discovery. In this chapter, we outline the methodological advances in FEP+, including the OPLS3 force fields, the REST2 (Replica Exchange with Solute Tempering) enhanced sampling, the incorporation of REST2 sampling with conventional FEP (Free Energy Perturbation) through FEP/REST, and the advanced simulation setup and data analysis. The validation of FEP+ method in retrospective studies and the prospective applications in drug discovery projects are also discussed. We then present the recent extension of FEP+ method to handle challenging perturbations, including core-hopping transformations, macrocycle modifications, and reversible covalent inhibitor optimization. The limitations and pitfalls of the current FEP+ methodology and the best practices in real applications are also examined.

Rigorous Free Energy Perturbation Approach to Estimating Relative Binding Affinities between Ligands with Multiple Protonation and Tautomeric States

Accurate prediction of ligand binding affinities is of key importance in small molecule lead optimization and a central task in computational medicinal chemistry. Over the years, advances in both computer hardware and computational methodologies have established free energy perturbation (FEP) methods as among the most reliable and rigorous approaches to compute protein-ligand binding free energies. However, accurate description of ionization and tautomerism of ligands is still a major challenge in structure-based prediction of binding affinities. Druglike molecules are often weak acid or bases with multiple accessible protonation and tautomeric states that can contribute significantly to the binding process. To address this issue, we introduce in this work the p K_a and tautomeric state correction approach. This approach is based on free energy perturbation formalism and provides a rigorous treatment of the ionization and tautomeric equilibria of ligands in solution and in the protein complexes. A series of Kinesin Spindle Protein (KSP) and Factor Xa inhibitor molecules were used as test cases. Our results demonstrate that the p K_a and tautomeric state correction approach is able to rigorously and accurately incorporate multiple protonation and tautomeric states in the binding affinity calculations.

Modeling the value of predictive affinity scoring in preclinical drug discovery

Drug discovery is widely recognized to be a difficult and costly activity in large part due to the challenge of identifying chemical matter which simultaneously optimizes multiple properties, one of which is affinity for the primary biological target. Further, many of these properties are difficult to predict ahead of expensive and time-consuming compound synthesis and experimental testing. Here we highlight recent work to develop compound affinity prediction models, and extensively investigate the value such models may provide to preclinical drug discovery. We demonstrate that the ability of these models to improve the overall probability of success is crucially dependent on the shape of the error distribution, not just the root-mean-square error. In particular, while scoring more molecule ideas generally improves the probability of project success when the error distribution is Gaussian, fat-tail distributions such as a Cauchy distribution, can lead to a situation where scoring more ideas actually decreases the overall probability of success.

Efficient sampling of puckering states of monosaccharides through replica exchange with solute tempering and bond softening

A molecular-level understanding of the structure, dynamics, and reactivity of carbohydrates is fundamental to the understanding of a range of key biological processes. The six-membered pyranose ring, a central component of biological monosaccharides and carbohydrates, has many different puckering conformations, and the conformational free energy landscape of these biologically important monosaccharides remains elusive. The puckering conformations of monosaccharides are separated by high energy barriers, which pose a great challenge for the complete sampling of these important conformations and accurate modeling of these systems. While metadynamics or umbrella sampling methods have been used to study the conformational space of monosaccharides, these methods might be difficult to generalize to other complex ring systems with more degrees of freedom. In this paper, we introduce a new enhanced sampling method for the rapid sampling over high energy barriers that combines our previously developed enhanced sampling method REST (replica exchange with solute tempering) with a bond softening (BOS) scheme that makes a chemical bond in the ring weaker as one ascends the replica ladder. We call this new method replica exchange with solute tempering and bond softening (REST/BOS). We demonstrate the superior sampling efficiency of REST/BOS over other commonly used enhanced sampling methods, including temperature replica exchange method and REST. The conformational free energy landscape of four biologically important monosaccharides, namely, α-glucose, β-glucose, β-mannose, and β-xylose, is studied using REST/BOS, and results are compared with previous experimental and theoretical studies.

High throughput evaluation of macrocyclization strategies for conformer stabilization

While macrocyclization of a linear compound to stabilize a known bioactive conformation can be a useful strategy to increase binding potency, the difficulty of macrocycle synthesis can limit the throughput of such strategies. Thus computational techniques may offer the higher throughput required to screen large numbers of compounds. Here we introduce a method for evaluating the propensity of a macrocyclic compound to adopt a conformation similar that of a known active linear compound in the binding site. This method can be used as a fast screening tool for prioritizing macrocycles by leveraging the assumption that the propensity for the known bioactive substructural conformation relates to the affinity. While this method cannot to identify new interactions not present in the known linear compound, it could quickly differentiate compounds where the three dimensional geometries imposed by the macrocyclization prevent adoption of conformations with the same contacts as the linear compound in their conserved region. Here we report the implementation of this method using an RMSD-based structural descriptor and a Boltzmann-weighted propensity calculation and apply it retrospectively to three macrocycle linker optimization design projects. We found the method performs well in terms of prioritizing more potent compounds.