Relative Binding Free Energy Calculations in Drug Discovery: Recent Advances and Practical Considerations

In Silico Predictions of Endocrine Disruptors Properties

Endocrine-disrupting chemicals (EDCs) are a broad class of molecules present in our environment that are suspected to cause adverse effects in the endocrine system by interfering with the synthesis, transport, degradation, or action of endogenous ligands. The characterization of the harmful interaction between environmental compounds and their potential cellular targets and the development of robust in vivo, in vitro, and in silico screening methods are important for assessment of the toxic potential of large numbers of chemicals. In this context, computer-aided technologies that will allow for activity prediction of endocrine disruptors and environmental risk assessments are being developed. These technologies must be able to cope with diverse data and connect chemistry at the atomic level with the biological activity at the cellular, organ, and organism levels. Quantitative structure-activity relationship methods became popular for toxicity issues. They correlate the chemical structure of compounds with biological activity through a number of molecular descriptors (e.g., molecular weight and parameters to account for hydrophobicity, topology, or electronic properties). Chemical structure analysis is a first step; however, modeling intermolecular interactions and cellular behavior will also be essential. The increasing number of three-dimensional crystal structures of EDCs' targets has provided a wealth of structural information that can be used to predict their interactions with EDCs using docking and scoring procedures. In the present review, we have described the various computer-assisted approaches that use ligands and targets properties to predict endocrine disruptor activities.

Optimal Measurement Network of Pairwise Differences

When both the difference between two quantities and their individual values can be measured or computationally predicted, multiple quantities can be determined from the measurements or predictions of select individual quantities and select pairwise differences. These measurements and predictions form a network connecting the quantities through their differences. Here, I analyze the optimization of such networks, where the trace (A-optimal), the largest eigenvalue (E-optimal), or the determinant (D-optimal) of the covariance matrix associated with the estimated quantities are minimized with respect to the allocation of the measurement (or computational) cost to different measurements (or predictions). My statistical analysis of the performance of such optimal measurement networks-based on large sets of simulated data-suggests that they substantially accelerate the determination of the quantities and that they may be useful in applications such as the computational prediction of binding free energies of candidate drug molecules.

Absolute Free Energy of Binding Calculations for Macrophage Migration Inhibitory Factor in Complex with a Druglike Inhibitor

Calculation of the absolute free energy of binding (ΔG_bind) for a complex in solution is challenging owing to the need for adequate configurational sampling and an accurate energetic description, typically with a force field (FF). In this study, Monte Carlo (MC) simulations with improved side-chain and backbone sampling are used to assess ΔG_bind for the complex of a druglike inhibitor (MIF180) with the protein macrophage migration inhibitory factor (MIF) using free energy perturbation (FEP) calculations. For comparison, molecular dynamics (MD) simulations were employed as an alternative sampling method for the same system. With the OPLS-AA/M FF and CM5 atomic charges for the inhibitor, the ΔG_bind results from the MC/FEP and MD/FEP simulations, -8.80 ± 0.74 and -8.46 ± 0.85 kcal/mol, agree well with each other and with the experimental value of -8.98 ± 0.28 kcal/mol. The convergence of the results and analysis of the trajectories indicate that sufficient sampling was achieved for both approaches. Repeating the MD/FEP calculations using current versions of the CHARMM and AMBER FFs led to a 6 kcal/mol range of computed ΔG_bind. These results show that calculation of accurate ΔG_bind for large ligands is both feasible and numerically equivalent, within error limits, using either methodology.

Exploring the Effectiveness of Binding Free Energy Calculations

Increasing the accuracy of the evaluation of ligand-binding energies is one of the most important tasks of current computational biology. Here we explore the accuracy of free energy perturbation (FEP) approaches by comparing the performance of a "regular" FEP method to the one using replica exchange to enhance the sampling on a well-defined benchmark. The examination was limited to the so-called alchemical perturbations which are restricted to a fragment of the drug, and therefore, the calculation is a relative one rather than the absolute binding energy of the drug. Overall, our calculations reach the 1 kcal/mol accuracy limit. It is also shown that the accurate prediction of the position of water molecules around the binding pocket is important for FEP calculations. Interestingly, the replica exchange method does not significantly improve the accuracy of binding energies, suggesting that we reach the limit where the force field quality is a critical factor for accurate calculations.

Target-Specific Prediction of Ligand Affinity with Structure-Based Interaction Fingerprints

Discovery and optimization of small molecule inhibitors as therapeutic drugs have immensely benefited from rational structure-based drug design. With recent advances in high-resolution structure determination, computational power, and machine learning methodology, it is becoming more tractable to elucidate the structural basis of drug potency. However, the applicability of machine learning models to drug design is limited by the interpretability of the resulting models in terms of feature importance. Here, we take advantage of the large number of available inhibitor-bound HIV-1 protease structures and associated potencies to evaluate inhibitor diversity and machine learning models to predict ligand affinity. First, using a hierarchical clustering approach, we grouped HIV-1 protease inhibitors and identified distinct core structures. Explicit features including protein-ligand interactions were extracted from high-resolution cocrystal structures as 3D-based fingerprints. We found that a gradient boosting machine learning model with this explicit feature attribution can predict binding affinity with high accuracy. Finally, Shapley values were derived to explain local feature importance. We found specific van der Waals (vdW) interactions of key protein residues are pivotal for the predicted potency. Protein-specific and interpretable prediction models can guide the optimization of many small molecule drugs for improved potency.

Development of a Robust Indirect Approach for MM → QM Free Energy Calculations That Combines Force-Matched Reference Potential and Bennett's Acceptance Ratio Methods

We use the PBE0/6-31G* density functional method to perform ab initio quantum mechanical/molecular mechanical (QM/MM) molecular dynamics (MD) simulations under periodic boundary conditions with rigorous electrostatics using the ambient potential composite Ewald method in order to test the convergence of MM → QM/MM free energy corrections for the prediction of 17 small-molecule solvation free energies and eight ligand binding free energies to T4 lysozyme. The "indirect" thermodynamic cycle for calculating free energies is used to explore whether a series of reference potentials improve the statistical quality of the predictions. Specifically, we construct a series of reference potentials that optimize a molecular mechanical (MM) force field's parameters to reproduce the ab initio QM/MM forces from a QM/MM simulation. The optimizations form a systematic progression of successively expanded parameters that include bond, angle, dihedral, and charge parameters. For each reference potential, we calculate benchmark quality reference values for the MM → QM/MM correction by performing the mixed MM and QM/MM Hamiltonians at 11 intermediate states, each for 200 ps. We then compare forward and reverse application of Zwanzig's relation, thermodynamic integration (TI), and Bennett's acceptance ratio (BAR) methods as a function of reference potential, simulation time, and the number of simulated intermediate states. We find that Zwanzig's equation is inadequate unless a large number of intermediate states are explicitly simulated. The TI and BAR mean signed errors are very small even when only the end-state simulations are considered, and the standard deviations of the TI and BAR errors are decreased by choosing a reference potential that optimizes the bond and angle parameters. We find a robust approach for the data sets of fairly rigid molecules considered here is to use bond + angle reference potential together with the end-state-only BAR analysis. This requires QM/MM simulations to be performed in order to generate reference data to parametrize the bond + angle reference potential, and then this same simulation serves a dual purpose as the full QM/MM end state. The convergence of the results with respect to time suggests that computational resources may be used more efficiently by running multiple simulations for no more than 50 ps, rather than running one long simulation.

Molecular Dynamics of Water Embedded Carbon Nanocones: Surface Waves Observation

We employed molecular dynamics simulations on the water solvation of conically shaped carbon nanoparticles. We explored the hydrophobic behaviour of the nanoparticles and investigated microscopically the cavitation of water in a conical confinement with different angles. We performed additional molecular dynamics simulations in which the carbon structures do not interact with water as if they were in vacuum. We detected a waving on the surface of the cones that resembles the shape agitations of artificial water channels and biological porins. The surface waves were induced by the pentagonal carbon rings (in an otherwise hexagonal network of carbon rings) concentrated near the apex of the cones. The waves were affected by the curvature gradients on the surface. They were almost undetected for the case of an armchair nanotube. Understanding such nanoscale phenomena is the key to better designed molecular models for membrane systems and nanodevices for energy applications and separation.

Overcoming Challenging Substituent Perturbations with Multisite λ-Dynamics: A Case Study Targeting β-Secretase 1

Alchemical free energy calculations have made a dramatic impact upon the field of structure-based drug design by allowing functional group modifications to be explored computationally prior to experimental synthesis and assay evaluation, thereby informing and directing synthetic strategies. In furthering the advancement of this area, a series of 21 β-secretase 1 (BACE1) inhibitors developed by Janssen Pharmaceuticals were examined to evaluate the ability to explore large substituent perturbations, some of which contain scaffold modifications, with multisite λ-dynamics (MSλD), an innovative alchemical free energy framework. Our findings indicate that MSλD is able to efficiently explore all structurally diverse ligand end-states simultaneously within a single MD simulation with a high degree of precision and with reduced computational costs compared to the widely used approach TI/MBAR. Furthermore, computational predictions were shown to be accurate to within 0.5-0.8 kcal/mol when CM1A partial atomic charges were combined with CHARMM or OPLS-AA-based force fields, demonstrating that MSλD is force field independent and a viable alternative to FEP or TI approaches for drug design.

Predicting Kinase Inhibitor Resistance: Physics-Based and Data-Driven Approaches

Resistance to small molecule drugs often emerges in cancer cells, viruses, and bacteria as a result of the evolutionary pressure exerted by the therapy. Protein mutations that directly impair drug binding are frequently involved in resistance, and the ability to anticipate these mutations would be beneficial in drug development and clinical practice. Here, we evaluate the ability of three distinct computational methods to predict ligand binding affinity changes upon protein mutation for the cancer target Abl kinase. These structure-based approaches rely on first-principle statistical mechanics, mixed physics- and knowledge-based potentials, and machine learning, and were able to estimate binding affinity changes and identify resistant mutations with remarkable accuracy. We expect that these complementary approaches will enable the routine prediction of resistance-causing mutations in a variety of other target proteins.

Robust Free Energy Perturbation Protocols for Creating Molecules in Solution

Accurate methods to estimate free energies play an important role for studying diverse condensed-phase problems in chemistry and biochemistry. The most common methods used in conjunction with molecular dynamics (MD) and Monte Carlo statistical mechanics (MC) simulations are free energy perturbation (FEP) and thermodynamic integration (TI). For common applications featuring the conversion of one molecule to another, simulations are run in stages or multiple "λ-windows" to promote convergence of the results. For computation of absolute free energies of solvation or binding, calculations are needed in which the solute is typically annihilated in the solvent and in the complex. The present work addresses identification of optimal protocols for such calculations, specifically, the creation/annihilation of organic molecules in aqueous solution. As is common practice, decoupling of the perturbations for electrostatic and Lennard-Jones interactions was performed. Consistent with earlier reports, FEP calculations for molecular creations are much more efficient, while annihilations require many more windows and may converge to incorrect values. Strikingly, we find that as few as four windows may be adequate for creation calculations for solutes ranging from argon to ethylbenzene. For a larger druglike molecule, MIF180, which contains 22 non-hydrogen atoms and three rotatable bonds, 10 creation windows are found to be adequate to yield the correct free energy of hydration. Convergence is impeded with procedures that use any sampling in the annihilation direction, and there is no need for postprocessing methods such as the Bennett acceptance ratio (BAR).

Computational Fluorine Scanning Using Free-Energy Perturbation

We present perturbative fluorine scanning, a computational fluorine scanning approach using free-energy perturbation. This method can be applied to molecular dynamics simulations of a single compound and make predictions for the best binders out of numerous fluorinated analogues. We tested the method on nine test systems: renin, DPP4, menin, P38, factor Xa, CDK2, AKT, JAK2, and androgen receptor. The predictions were in excellent agreement with more rigorous alchemical free-energy calculations and in good agreement with experimental data for most of the test systems. However, the agreement with experiment was very poor in some of the test systems, and this highlights the need for improved force fields in addition to accurate treatment of tautomeric and protonation states. The method is of particular interest due to the wide use of fluorine in medicinal chemistry to improve binding affinity and ADME properties. The promising results on this test case suggest that perturbative fluorine scanning will be a useful addition to the available arsenal of free-energy methods.

Error Tolerance of Machine Learning Algorithms across Contemporary Biological Targets

Machine learning continues to make strident advances in the prediction of desired properties concerning drug development. Problematically, the efficacy of machine learning in these arenas is reliant upon highly accurate and abundant data. These two limitations, high accuracy and abundance, are often taken together; however, insight into the dataset accuracy limitation of contemporary machine learning algorithms may yield insight into whether non-bench experimental sources of data may be used to generate useful machine learning models where there is a paucity of experimental data. We took highly accurate data across six kinase types, one GPCR, one polymerase, a human protease, and HIV protease, and intentionally introduced error at varying population proportions in the datasets for each target. With the generated error in the data, we explored how the retrospective accuracy of a Naïve Bayes Network, a Random Forest Model, and a Probabilistic Neural Network model decayed as a function of error. Additionally, we explored the ability of a training dataset with an error profile resembling that produced by the Free Energy Perturbation method (FEP+) to generate machine learning models with useful retrospective capabilities. The categorical error tolerance was quite high for a Naïve Bayes Network algorithm averaging 39% error in the training set required to lose predictivity on the test set. Additionally, a Random Forest tolerated a significant degree of categorical error introduced into the training set with an average error of 29% required to lose predictivity. However, we found the Probabilistic Neural Network algorithm did not tolerate as much categorical error requiring an average of 20% error to lose predictivity. Finally, we found that a Naïve Bayes Network and a Random Forest could both use datasets with an error profile resembling that of FEP+. This work demonstrates that computational methods of known error distribution like FEP+ may be useful in generating machine learning models not based on extensive and expensive in vitro-generated datasets.

Template-Based Method for Conformation Generation and Scoring for Congeneric Series of Ligands

Physics-based prediction of protein-ligand binding affinities for a congeneric series of ligands in lead optimization requires their geometries as a first step. In this paper, we report a method that uses the 3D conformation of a lead compound in complex with a protein as a template to generate conformations of a series of related analog compounds. The method uses the Maximal Common Substructure (MCS) computed between lead and analog ligands to assign coordinates for the atoms shared between the ligands. For the differing atoms, a conformation generation procedure is implemented that results in a diversity of conformations. The generated conformations are sorted using a score based on the Molecular Mechanics and Generalized Born with Solvent Accessible Surface Area contribution (MM-GBSA) method. The accuracy of the generated conformations is tested retrospectively using a cross-validation approach applied to four data sets obtained from the Drug Design Data Resource (D3R) by measuring the RMSD of the top scored conformation with respect to the crystallographic pose. The scoring ability of the method is independently assessed using data for the same protein targets to test the rank ordering ability and separating active and inactive ligands. We tested the effect of protein flexibility during structural optimization and scoring approaches with and without strain energies. Retrospective validation on data sets comprising 4 targets shows that the method outperforms random selection for all targets and outperforms a molecular weight-based null model in 3 out of 4 targets in separating active and inactive compounds. Therefore, the presented method is expected to be of utility in lead optimization for rapidly screening analog ligands and generating initial conformations for use in more detailed physics-based binding affinity prediction methods.

Three dimensional structure prediction of panomycocin, a novel Exo-β-1,3-glucanase isolated from Wickerhamomyces anomalus NCYC 434 and the computational site-directed mutagenesis studies to enhance its thermal stability for therapeutic applications

Panomycocin is a naturally produced potent antimycotic/antifungal protein secreted by the yeast Wickerhamomyces anomalus NCYC 434 with an exo-β-1,3-glucanase activity. In this study the three dimensional structure of panomycocin was predicted and the computational site-directed mutagenesis was performed to enhance its thermal stability in liquid formulations over the body temperature for topical therapeutic applications. Homology modeling was performed with MODELLER and I-TASSER. Among the generated models, the model with the lowest energy and DOPE score was selected for further loop modeling. The loop model was optimized and the reliability of the model was confirmed with ERRAT, Verify 3D and Ramachandran plot values. Enhancement of the thermal stability of the model was done using contemporary servers and programs such as SPDBViewer, CNA, I-Mutant2.0, Eris, AUTO-MUTE and MUpro. In the region outside the binding site of the model Leu52 Arg, Phe223Arg and Gly254Arg were found to be the best thermostabilizing mutations with 6.26 K, 6.26 K and 8.27 K increases, respectively. In the binding site Glu186Arg was found to be the best thermostabilizer mutation with a 9.58 K temperature increase. The results obtained in this study led us to design a mutant panomycocin that can be used as a novel antimycotic/antifungal drug in a liquid formulation for topical applications over the normal body temperature.

Molecular Modeling-Guided Design of Phospholipid-Based Prodrugs

The lipidic prodrug approach is an emerging field for improving a number of biopharmaceutical and drug delivery aspects. Owing to their structure and nature, phospholipid (PL)-based prodrugs may join endogenous lipid processing pathways, and hence significantly improve the pharmacokinetics and/or bioavailability of the drug. Additional advantages of this approach include drug targeting by enzyme-triggered drug release, blood–brain barrier permeability, lymphatic targeting, overcoming drug resistance, or enabling appropriate formulation. The PL-prodrug design includes various structural modalities-different conjugation strategies and/or the use of linkers between the PL and the drug moiety, which considerably influence the prodrug characteristics and the consequent effects. In this article, we describe how molecular modeling can guide the structural design of PL-based prodrugs. Computational simulations can predict the extent of phospholipase A₂ (PLA₂)-mediated activation, and facilitate prodrug development. Several computational methods have been used to facilitate the design of the pro-drugs, which will be reviewed here, including molecular docking, the free energy perturbation method, molecular dynamics simulations, and free density functional theory. Altogether, the studies described in this article indicate that computational simulation-guided PL-based prodrug molecular design correlates well with the experimental results, allowing for more mechanistic and less empirical development. In the future, the use of molecular modeling techniques to predict the activity of PL-prodrugs should be used earlier in the development process.

An Overview of Molecular Modeling for Drug Discovery with Specific Illustrative Examples of Applications

In this paper we review the current status of high-performance computing applications in the general area of drug discovery. We provide an introduction to the methodologies applied at atomic and molecular scales, followed by three specific examples of implementation of these tools. The first example describes in silico modeling of the adsorption of small molecules to organic and inorganic surfaces, which may be applied to drug delivery issues. The second example involves DNA translocation through nanopores with major significance to DNA sequencing efforts. The final example offers an overview of computer-aided drug design, with some illustrative examples of its usefulness.

QligFEP: an automated workflow for small molecule free energy calculations in Q

The process of ligand binding to a biological target can be represented as the equilibrium between the relevant solvated and bound states of the ligand. This which is the basis of structure-based, rigorous methods such as the estimation of relative binding affinities by free energy perturbation (FEP). Despite the growing capacity of computing power and the development of more accurate force fields, a high throughput application of FEP is currently hampered due to the need, in the current schemes, of an expert user definition of the "alchemical" transformations between molecules in the series explored. Here, we present QligFEP, a solution to this problem using an automated workflow for FEP calculations based on a dual topology approach. In this scheme, the starting poses of each of the two ligands, for which the relative affinity is to be calculated, are explicitly present in the MD simulations associated with the (dual topology) FEP transformation, making the perturbation pathway between the two ligands univocal. We show that this generalized method can be applied to accurately estimate solvation free energies for amino acid sidechain mimics, as well as the binding affinity shifts due to the chemical changes typical of lead optimization processes. This is illustrated in a number of protein systems extracted from other FEP studies in the literature: inhibitors of CDK2 kinase and a series of A2A adenosine G protein-coupled receptor antagonists, where the results obtained with QligFEP are in excellent agreement with experimental data. In addition, our protocol allows for scaffold hopping perturbations to identify the binding affinities between different core scaffolds, which we illustrate with a series of Chk1 kinase inhibitors. QligFEP is implemented in the open-source MD package Q, and works with the most common family of force fields: OPLS, CHARMM and AMBER.

Effect of set up protocols on the accuracy of alchemical free energy calculation over a set of ACK1 inhibitors

Hit-to-lead virtual screening frequently relies on a cascade of computational methods that starts with rapid calculations applied to a large number of compounds and ends with more expensive computations restricted to a subset of compounds that passed initial filters. This work focuses on set up protocols for alchemical free energy (AFE) scoring in the context of a Docking–MM/PBSA–AFE cascade. A dataset of 15 congeneric inhibitors of the ACK1 protein was used to evaluate the performance of AFE set up protocols that varied in the steps taken to prepare input files (using previously docked and best scored poses, manual selection of poses, manual placement of binding site water molecules). The main finding is that use of knowledge derived from X-ray structures to model binding modes, together with the manual placement of a bridging water molecule, improves the R² from 0.45 ± 0.06 to 0.76 ± 0.02 and decreases the mean unsigned error from 2.11 ± 0.08 to 1.24 ± 0.04 kcal mol^-1. By contrast a brute force automated protocol that increased the sampling time ten-fold lead to little improvements in accuracy. Besides, it is shown that for the present dataset hysteresis can be used to flag poses that need further attention even without prior knowledge of experimental binding affinities.

Assessing the performance of three resveratrol in binding with SIRT1 by molecular dynamics simulation and MM/GBSA methods: the weakest binding of resveratrol 3 to SIRT1 triggers a possibility of dissociation from its binding site

SIRT1 is an NAD+-dependent deacetylase, whose activators have potential therapeutic applications in the age-related, metabolic, neurode-generative and cardiovascular diseases. Resveratrol (RSV) has been regarded as potent SIRT1 activator in the treatment of atherosclerosis and cardiovascular diseases. Moreover, the previous crystal structure of SIRT1 complex with RSV indicated that three RSV participated in the SIRT1 activation, namely RSV1 and RSV2 strengthened the interaction of SIRT1 to its substrate peptide and promoted the stimulation of SIRT1 activity, while RSV3 exhibited an auxiliary function. Whereas the molecular mechanism of the insignificant role of RSV3 has not been reported. Moreover, the dynamics properties of three RSV in the binding site of SIRT1 are still unknown. In this study, to evaluate the role of RSV3 systematically, five SIRT1_RSV system (SIRT1_RSV^1-3, SIRT1_RSV^1-2, SIRT1_RSV^{1, 3}, SIRT1_RSV^2-3 and SIRT1³) were constructed and subjected to 20 ns molecular dynamics simulation. Our simulations showed that the binding of RSV1 was extremely stable in any system and might be independent to the existence of RSV2 and RSV3. The stability of RSV2 was slightly weaker than that of RSV1. More interestingly, RSV3 could completely dissociate from its binding site in SIRT1_RSV^1-3 system at 7 ns, which might be an important factor for its auxiliary role in activating SIRT1. Actually, some preliminary signs of RSV3 dissociation already appeared around 1 ns. In detail, the innate weakness between RSV3 and Lys444 would trigger frequent formation and breakage of the RSV3-Lys444 hydrogen bond. Residue Asp298 preferred to form hydrogen bond with RSV2 rather than RSV3 once the RSV3-Lys444 hydrogen bond became stronger. Without the RSV3-Asp298 interaction, the remaining RSV3-Lys444 and RSV3-Asp292 hydrogen bonds hardly preserved the RSV3's binding, which resulted in the dissociation of RSV3. Our work elucidated the inequable contributions of three RSV on their binding interaction with SIRT1 and indicated that the binding sites of RSV1 and RSV2 should be enough to design more potent SIRT1 activator.

Absolute Binding Free Energy Calculation and Design of a Subnanomolar Inhibitor of Phosphodiesterase-10

Accurate prediction of absolute protein-ligand binding free energy could considerably enhance the success rate of structure-based drug design but is extremely challenging and time-consuming. Free energy perturbation (FEP) has been proven reliable but is limited to prediction of relative binding free energies of similar ligands (with only minor structural differences) in binding with a same drug target in practical drug design applications. Herein, a Gaussian algorithm-enhanced FEP (GA-FEP) protocol has been developed to enhance the FEP simulation performance, enabling to efficiently carry out the FEP simulations on vanishing the whole ligand and, thus, predict the absolute binding free energies (ABFEs). Using the GA-FEP protocol, the FEP simulations for the ABFE calculation (denoted as GA-FEP/ABFE) can achieve a satisfactory accuracy for both structurally similar and diverse ligands in a dataset of more than 100 receptor-ligand systems. Further, our GA-FEP/ABFE-guided lead optimization against phosphodiesterase-10 led to the discovery of a subnanomolar inhibitor (IC₅₀ = 0.87 nM, ∼2000-fold improvement in potency) with cocrystal confirmation.

The nature of ligand efficiency

Ligand efficiency is a widely used design parameter in drug discovery. It is calculated by scaling affinity by molecular size and has a nontrivial dependency on the concentration unit used to express affinity that stems from the inability of the logarithm function to take dimensioned arguments. Consequently, perception of efficiency varies with the choice of concentration unit and it is argued that the ligand efficiency metric is not physically meaningful nor should it be considered to be a metric. The dependence of ligand efficiency on the concentration unit can be eliminated by defining efficiency in terms of sensitivity of affinity to molecular size and this is illustrated with reference to fragment-to-lead optimizations. Group efficiency and fit quality are also examined in detail from a physicochemical perspective. The importance of examining relationships between affinity and molecular size directly is stressed throughout this study and an alternative to ligand efficiency for normalization of affinity with respect to molecular size is presented.

A computationally designed binding mode flip leads to a novel class of potent tri-vector cyclophilin inhibitors

Cyclophilins (Cyps) are a major family of drug targets that are challenging to prosecute with small molecules because the shallow nature and high degree of conservation of the active site across human isoforms offers limited opportunities for potent and selective inhibition. Herein a computational approach based on molecular dynamics simulations and free energy calculations was combined with biophysical assays and X-ray crystallography to explore a flip in the binding mode of a reported urea-based Cyp inhibitor. This approach enabled access to a distal pocket that is poorly conserved among key Cyp isoforms, and led to the discovery of a new family of sub-micromolar cell-active inhibitors that offer unprecedented opportunities for the development of next-generation drug therapies based on Cyp inhibition. The computational approach is applicable to a broad range of organic functional groups and could prove widely enabling in molecular design.

Characterizing early drug resistance-related events using geometric ensembles from HIV protease dynamics

The use of antiretrovirals (ARVs) has drastically improved the life quality and expectancy of HIV patients since their introduction in health care. Several millions are still afflicted worldwide by HIV and ARV resistance is a constant concern for both healthcare practitioners and patients, as while treatment options are finite, the virus constantly adapts via complex mutation patterns to select for resistant strains under the pressure of drug treatment. The HIV protease is a crucial enzyme for viral maturation and has been a game changing drug target since the first application. Due to similarities in protease inhibitor designs, drug cross-resistance is not uncommon across ARVs of the same class. It is known that resistance against protease inhibitors is associated with a wider active site, but results from our large scale molecular dynamics simulations combined with statistical tests and network analysis further show, for the first time, that there are regions of local expansions and compactions associated with high levels of resistance conserved across eight different protease inhibitors visible in their complexed form within closed receptor conformations. The observed conserved expansion sites may provide an alternative drug-targeting site. Further, the method developed here is novel, supplementary to methods of variation analysis at sequence level, and should be applicable in analysing the structural consequences of mutations in other contexts using molecular ensembles.

Systematic Mutagenesis of Oncocin Reveals Enhanced Activity and Insights into the Mechanisms of Antimicrobial Activity

Oncocin is a proline-rich antimicrobial peptide that inhibits protein synthesis by binding to the bacterial ribosome. In this work, the antimicrobial activity of oncocin was improved by systematic peptide mutagenesis and activity evaluation. We found that a pair of cationic substitutions (P4K and L7K/R) improves the activity by 2-4 fold (p<0.05) against multiple Gram-negative bacteria. An in vitro transcription / translation assay indicated that the increased activity was not because of stronger ribosome binding. Rather a cellular internalization assay revealed a higher internalization rate for the optimized analogs thereby suggesting a mechanism to increase potency. In addition, we found that the optimized peptides' benefit is dependent upon nutrient-depleted media conditions. The molecular design and characterization strategies have broad potential for development of antimicrobial peptides.

A survey of the role of nitrile groups in protein–ligand interactions

In classical medicinal chemistry, nitrile groups were commonly considered as bioisosteres of carbonyl, hydroxyl and carboxyl groups, as well as halogen atoms. However, there is a lack of in-depth understanding about the structural and energetic characteristics of nitrile groups in protein–ligand interactions. Here, we have surveyed the Protein Data Bank and ChEMBL databases with the goal of characterizing such protein–ligand interactions for nitrile-containing compounds. We discuss the versatile roles of nitrile groups in improving binding affinities, and give special attention to examples of displacing and mimicking binding-site waters by nitrile groups. We expect that this review article will further inspire medicinal chemists to exploit nitrile groups rationally in structure-based drug design.

Ligand Binding Free Energies with Adaptive Water Networks: Two-Dimensional Grand Canonical Alchemical Perturbations

Computational methods to calculate ligand binding affinities are increasing in popularity, due to improvements in simulation algorithms, computational resources, and easy-to-use software. However, issues can arise in relative ligand binding free energy simulations if the ligands considered have different active site water networks, as simulations are typically performed with a predetermined number of water molecules (fixed N ensembles) in preassigned locations. If an alchemical perturbation is attempted where the change should result in a different active site water network, the water molecules may not be able to adapt appropriately within the time scales of the simulations—particularly if the active site is occluded. By combining the grand canonical ensemble (μVT) to sample active site water molecules, with conventional alchemical free energy methods, the water network is able to dynamically adapt to the changing ligand. We refer to this approach as grand canonical alchemical perturbation (GCAP). In this work we demonstrate GCAP for two systems; Scytalone Dehydratase (SD) and Adenosine A_2A receptor. For both systems, GCAP is shown to perform well at reproducing experimental binding affinities. Calculating the relative binding affinities with a naïve, conventional attempt to solvate the active site illustrates how poor results can be if proper consideration of water molecules in occluded pockets is neglected. GCAP results are shown to be consistent with time-consuming double decoupling simulations. In addition, by obtaining the free energy surface for ligand perturbations, as a function of both the free energy coupling parameter and water chemical potential, it is possible to directly deconvolute the binding energetics in terms of protein–ligand direct interactions and protein binding site hydration.

From bench to bedside, via desktop. Recent advances in the application of cutting-edge in silico tools in the research of drugs targeting bromodomain modules

The discipline of drug discovery has greatly benefited by computational tools and in silico algorithms aiming at rationalization of many related processes, from the stage of early hit identification to the preclinical phases of drug candidate validation. The various methodologies referred to as molecular modeling tools span a broad spectrum of applications, from straightforward approaches such as virtual screening of compound libraries to more advanced techniques involving the precise estimation of free energy upon binding of the candidate drug to its macromolecular target. To this end, we report an overview of specific studies where implementation of such sophisticated modeling algorithms has shown to be indispensable for addressing challenging systems and biological questions otherwise difficult to answer. We focus our attention on the emerging field of bromodomain inhibitors. Bromodomains are small modules involved in epigenetic signaling and currently comprise high-priority targets for developing both drug candidates and chemical probes for basic biomedical research. We attempt a critical presentation of selected cases utilizing cutting-edge in silico methodologies, with our main emphasis being on absolute or relative free energy simulations, on implementation of quantum-mechanics level calculations and on characterization of solvent thermodynamics. We discuss the advantages and strengths as well as the drawbacks and weaknesses of computational tools utilized in those works and we attempt to comment on specific issues related to their integration into the regular medicinal chemistry practice. Our conclusion is that while such methods indeed represent highly promising resources for further advancing the discipline, their application is not always trivial.

Reproducibility of Free Energy Calculations across Different Molecular Simulation Software Packages

Alchemical free energy calculations are an increasingly important modern simulation technique to calculate free energy changes on binding or solvation. Contemporary molecular simulation software such as AMBER, CHARMM, GROMACS, and SOMD include support for the method. Implementation details vary among those codes, but users expect reliability and reproducibility, i.e., for a given molecular model and set of force field parameters, comparable free energy differences should be obtained within statistical bounds regardless of the code used. Relative alchemical free energy (RAFE) simulation is increasingly used to support molecule discovery projects, yet the reproducibility of the methodology has been less well tested than its absolute counterpart. Here we present RAFE calculations of hydration free energies for a set of small organic molecules and demonstrate that free energies can be reproduced to within about 0.2 kcal/mol with the aforementioned codes. Absolute alchemical free energy simulations have been carried out as a reference. Achieving this level of reproducibility requires considerable attention to detail and package-specific simulation protocols, and no universally applicable protocol emerges. The benchmarks and protocols reported here should be useful for the community to validate new and future versions of software for free energy calculations.

Overview of the SAMPL6 host-guest binding affinity prediction challenge

Accurately predicting the binding affinities of small organic molecules to biological macromolecules can greatly accelerate drug discovery by reducing the number of compounds that must be synthesized to realize desired potency and selectivity goals. Unfortunately, the process of assessing the accuracy of current computational approaches to affinity prediction against binding data to biological macromolecules is frustrated by several challenges, such as slow conformational dynamics, multiple titratable groups, and the lack of high-quality blinded datasets. Over the last several SAMPL blind challenge exercises, host-guest systems have emerged as a practical and effective way to circumvent these challenges in assessing the predictive performance of current-generation quantitative modeling tools, while still providing systems capable of possessing tight binding affinities. Here, we present an overview of the SAMPL6 host-guest binding affinity prediction challenge, which featured three supramolecular hosts: octa-acid (OA), the closely related tetra-endo-methyl-octa-acid (TEMOA), and cucurbit[8]uril (CB8), along with 21 small organic guest molecules. A total of 119 entries were received from ten participating groups employing a variety of methods that spanned from electronic structure and movable type calculations in implicit solvent to alchemical and potential of mean force strategies using empirical force fields with explicit solvent models. While empirical models tended to obtain better performance than first-principle methods, it was not possible to identify a single approach that consistently provided superior results across all host-guest systems and statistical metrics. Moreover, the accuracy of the methodologies generally displayed a substantial dependence on the system considered, emphasizing the need for host diversity in blind evaluations. Several entries exploited previous experimental measurements of similar host-guest systems in an effort to improve their physical-based predictions via some manner of rudimentary machine learning; while this strategy succeeded in reducing systematic errors, it did not correspond to an improvement in statistical correlation. Comparison to previous rounds of the host-guest binding free energy challenge highlights an overall improvement in the correlation obtained by the affinity predictions for OA and TEMOA systems, but a surprising lack of improvement regarding root mean square error over the past several challenge rounds. The data suggests that further refinement of force field parameters, as well as improved treatment of chemical effects (e.g., buffer salt conditions, protonation states), may be required to further enhance predictive accuracy.

Estimating Atomic Contributions to Hydration and Binding Using Free Energy Perturbation

We present a general method called atom-wise free energy perturbation (AFEP), which extends a conventional molecular dynamics free energy perturbation (FEP) simulation to give the contribution to a free energy change from each atom. AFEP is derived from an expansion of the Zwanzig equation used in the exponential averaging method by defining that the system total energy can be partitioned into contributions from each atom. A partitioning method is assumed and used to group terms in the expansion to correspond to individual atoms. AFEP is applied to six example free energy changes to demonstrate the method. Firstly, the hydration free energies of methane, methanol, methylamine, methanethiol, and caffeine in water. AFEP highlights the atoms in the molecules that interact favorably or unfavorably with water. Finally AFEP is applied to the binding free energy of human immunodeficiency virus type 1 protease to lopinavir, and AFEP reveals the contribution of each atom to the binding free energy, indicating candidate areas of the molecule to improve to produce a more strongly binding inhibitor. FEP gives a single value for the free energy change and is already a very useful method. AFEP gives a free energy change for each "part" of the system being simulated, where part can mean individual atoms, chemical groups, amino acids, or larger partitions depending on what the user is trying to measure. This method should have various applications in molecular dynamics studies of physical, chemical, or biochemical phenomena, specifically in the field of computational drug discovery.

Predicting Binding Free Energies of PDE2 Inhibitors. The Difficulties of Protein Conformation

A congeneric series of 21 phosphodiesterase 2 (PDE2) inhibitors are reported. Crystal structures show how the molecules can occupy a 'top-pocket' of the active site. Molecules with small substituents do not enter the pocket, a critical leucine (Leu770) is closed and water molecules are present. Large substituents enter the pocket, opening the Leu770 conformation and displacing the waters. We also report an X-ray structure revealing a new conformation of the PDE2 active site domain. The relative binding affinities of these compounds were studied with free energy perturbation (FEP) methods and it represents an attractive real-world test case. In general, the calculations could predict the energy of small-to-small, or large-to-large molecule perturbations. However, accurately capturing the transition from small-to-large proved challenging. Only when using alternative protein conformations did results improve. The new X-ray structure, along with a modelled dimer, conferred stability to the catalytic domain during the FEP molecular dynamics (MD) simulations, increasing the convergence and thereby improving the prediction of ΔΔG of binding for some small-to-large transitions. In summary, we found the most significant improvement in results when using different protein structures, and this data set is useful for future free energy validation studies.