Performance of Relative Binding Free Energy Calculations on an Automatically Generated Dataset of Halogen-Deshalogen Matched Molecular Pairs.

In this study, we generated a matched molecular pair dataset of halogen/deshalogen compounds with reliable binding affinity data and structural binding mode information from public databases. The workflow includes automated system preparation and setup of free energy perturbation relative binding free energy calculations. We demonstrate the suitability of these datasets to investigate the performance of molecular mechanics force fields and molecular simulation algorithms for the purpose of in silico affinity predictions in lead optimization. Our datasets of a total of 115 matched molecular pairs show highly accurate binding free energy predictions with an average error of <1 kcal/mol despite the semi-automated calculation scheme. We quantify the accuracy of the optimized potential for liquid simulations (OPLS) force field to predict the effect of halogen addition to compounds, a commonly employed chemical modification in the design of drug-like molecules.

Large-Scale Assessment of Binding Free Energy Calculations in Active Drug Discovery Projects.

Accurate ranking of compounds with regards to their binding affinity to a protein using computational methods is of great interest to pharmaceutical research. Physics-based free energy calculations are regarded as the most rigorous way to estimate binding affinity. In recent years, many retrospective studies carried out both in academia and industry have demonstrated its potential. Here, we present the results of large-scale prospective application of the FEP+ method in active drug discovery projects in an industry setting at Merck KGaA, Darmstadt, Germany. We compare these prospective data to results obtained on a new diverse, public benchmark of eight pharmaceutically relevant targets. Our results offer insights into the challenges faced when using free energy calculations in real-life drug discovery projects and identify limitations that could be tackled by future method development. The new public data set we provide to the community can support further method development and comparative benchmarking of free energy calculations.

Large-scale in vitro functional testing and novel variant scoring via protein modeling provide insights into alkaline phosphatase activity in hypophosphatasia.

Hypophosphatasia (HPP) is a rare metabolic disorder characterized by low tissue-nonspecific alkaline phosphatase (TNSALP) typically caused by ALPL gene mutations. HPP is heterogeneous, with clinical presentation correlating with residual TNSALP activity and/or dominant-negative effects (DNE). We measured residual activity and DNE for 155 ALPL variants by transient transfection and TNSALP enzymatic activity measurement. Ninety variants showed low residual activity and 24 showed DNE. These results encompass all missense variants with carrier frequencies above 1/25,000 from the Genome Aggregation Database. We used resulting data as a reference to develop a new computational algorithm that scores ALPL missense variants and predicts high/low TNSALP enzymatic activity. Our approach measures the effects of amino acid changes on TNSALP dimer stability with a physics-based implicit solvent energy model. We predict mutation deleteriousness with high specificity, achieving a true-positive rate of 0.63 with false-positive rate of 0, with an area under receiver operating curve (AUC) of 0.9, better than all in silico predictors tested. Combining this algorithm with other in silico approaches can further increase performance, reaching an AUC of 0.94. This study expands our understanding of HPP heterogeneity and genotype/phenotype relationships with the aim of improving clinical ALPL variant interpretation.

Best of Two Worlds? How MD Simulations of Amphiphilic Helical Peptides in Membranes Can Complement Data from Oriented Solid-State NMR.

The membrane alignment of helical amphiphilic peptides in oriented phospholipid bilayers can be obtained as ensemble and time averages from solid state 2H NMR by fitting the quadrupolar splittings to ideal α-helices. At the same time, molecular dynamics (MD) simulations can provide atomistic insight into peptide-membrane systems. Here, we evaluate the potential of MD simulations to complement the experimental NMR data that is available on three exemplary systems: the natural antimicrobial peptide PGLa and the two designer-made peptides MSI-103 and KIA14, whose sequences were derived from PGLa. Each peptide was simulated for 1 µs in a DMPC lipid bilayer. We calculated from the MD simulations the local angles which define the side chain geometry with respect to the peptide helix. The peptide orientation was then calculated (i) directly from the simulation, (ii) from back-calculated MD-derived NMR splittings, and (iii) from experimental 2H NMR splittings. Our findings are that (1) the membrane orientation and secondary structure of the peptides found in the NMR analysis are generally well reproduced by the simulations; (2) the geometry of the side chains with respect to the helix backbone can deviate significantly from the ideal structure depending on the specific residue, but on average all side chains have the same orientation; and (3) for all of our peptides, the azimuthal rotation angle found from the MD-derived splittings is about 15° smaller than the experimental value.

Thermodynamic integration network approach to ion transport through protein channels: Perspectives and limits.

We present a molecular dynamics simulation study of alkali metal cation transport through the double-helical and the head-to-head conformers of the gramicidin ion channel. Our approach is based on a thermodynamic integration network, which consists of a sequence of transport reactions, absolute free energies of solvation and cycles of alchemical transmutations of the ions. In this manner, we can reliably estimate free energies and their statistical errors via a least-squares method without imposing external forces on the system. Within the double helical channel, we find a free energy surface typical for hopping transport between isoenergetic sites of ion localization, separated by comparatively large activation barriers. For fast transport through the head-to-head conformation, the thermodynamic network scheme starts to break down. © 2018 Wiley Periodicals, Inc.

Predicting resistance of clinical Abl mutations to targeted kinase inhibitors using alchemical free-energy calculations.

The therapeutic effect of targeted kinase inhibitors can be significantly reduced by intrinsic or acquired resistance mutations that modulate the affinity of the drug for the kinase. In cancer, the majority of missense mutations are rare, making it difficult to predict their impact on inhibitor affinity. This complicates the practice of precision medicine, pairing of patients with clinical trials, and development of next-generation inhibitors. Here, we examine the potential for alchemical free-energy calculations to predict how kinase mutations modulate inhibitor affinities to Abl, a major target in chronic myelogenous leukemia (CML). We find these calculations can achieve useful accuracy in predicting resistance for a set of eight FDA-approved kinase inhibitors across 144 clinically-identified point mutations, achieving a root mean square error in binding free energy changes of 1.10.91.3 kcal/mol (95% confidence interval) and correctly classifying mutations as resistant or susceptible with 888293% accuracy. Since these calculations are fast on modern GPUs, this benchmark establishes the potential for physical modeling to collaboratively support the rapid assessment and anticipation of the potential for patient mutations to affect drug potency in clinical applications.

Validation of tautomeric and protomeric binding modes by free energy calculations. A case study for the structure based optimization of D-amino acid oxidase inhibitors.

Optimization of fragment size D-amino acid oxidase (DAAO) inhibitors was investigated using a combination of computational and experimental methods. Retrospective free energy perturbation (FEP) calculations were performed for benzo[d]isoxazole derivatives, a series of known inhibitors with two potential binding modes derived from X-ray structures of other DAAO inhibitors. The good agreement between experimental and computed binding free energies in only one of the hypothesized binding modes strongly support this bioactive conformation. Then, a series of 1-H-indazol-3-ol derivatives formerly not described as DAAO inhibitors was investigated. Binding geometries could be reliably identified by structural similarity to benzo[d]isoxazole and other well characterized series and FEP calculations were performed for several tautomers of the deprotonated and protonated compounds since all these forms are potentially present owing to the experimental pKa values of representative compounds in the series. Deprotonated compounds are proposed to be the most important bound species owing to the significantly better agreement between their calculated and measured affinities compared to the protonated forms. FEP calculations were also used for the prediction of the affinities of compounds not previously tested as DAAO inhibitors and for a comparative structure-activity relationship study of the benzo[d]isoxazole and indazole series. Selected indazole derivatives were synthesized and their measured binding affinity towards DAAO was in good agreement with FEP predictions.

Prediction of Accurate Binding Modes Using Combination of Classical and Accelerated Molecular Dynamics and Free-Energy Perturbation Calculations: An Application to Toxicity Studies.

Estimating the correct binding modes of ligands in protein-ligand complexes is crucial not only in the drug discovery process but also for elucidating potential toxicity mechanisms. In the current paper, we propose a computational modeling workflow using the combination of docking, classical molecular dynamics (cMD), accelerated molecular dynamics (aMD) and free-energy perturbation (FEP+ protocol) for identification of possible ligand binding modes. It was applied for investigation of selected perfluorocarboxyl acids (PFCAs) in the PPARγ nuclear receptor. Although both regular and induced fit docking failed to reproduce the experimentally determined binding mode of the ligands when docked into a non-native X-ray structure, cMD and aMD simulations successfully identified the most probable binding conformations. Moreover, multiple binding modes were identified for all of these compounds and the shorter-chain PFCAs continuously moved between a few energetically favorable binding conformations. On the basis of MD predictions of binding conformations, we applied the default and also redesigned FEP+ sampling protocols, which accurately reproduced experimental differences in the binding energies. Thus, the preliminary MD simulations can also provide helpful information about correct setup of the FEP+ calculations. These results show that the PFCA binding modes were accurately predicted and that the FEP+ protocol can be used to estimate free energies of binding of flexible ligands that are not typical druglike compounds. Our in silico workflow revealed the specific ligand-residue interactions within the ligand binding domain and the main characteristics of the PFCAs, and it was concluded that these compounds are week PPARγ partial agonists. This work also suggests a common pipeline for identification of ligand binding modes, ligand-protein dynamics description, and relative free-energy calculations.

A Water-Bridged Cysteine-Cysteine Redox Regulation Mechanism in Bacterial Protein Tyrosine Phosphatases.

The emergence of multidrug-resistant Mycobacterium tuberculosis (Mtb) strains highlights the need to develop more efficacious and potent drugs. However, this goal is dependent on a comprehensive understanding of Mtb virulence protein effectors at the molecular level. Here, we used a post-expression cysteine (Cys)-to-dehydrolanine (Dha) chemical editing strategy to identify a water-mediated motif that modulates accessibility of the protein tyrosine phosphatase A (PtpA) catalytic pocket. Importantly, this water-mediated Cys-Cys non-covalent motif is also present in the phosphatase SptpA from Staphylococcus aureus, which suggests a potentially preserved structural feature among bacterial tyrosine phosphatases. The identification of this structural water provides insight into the known resistance of Mtb PtpA to the oxidative conditions that prevail within an infected host macrophage. This strategy could be applied to extend the understanding of the dynamics and function(s) of proteins in their native state and ultimately aid in the design of small-molecule modulators.

Automated Transition State Search and Its Application to Diverse Types of Organic Reactions.

Transition state search is at the center of multiple types of computational chemical predictions related to mechanistic investigations, reactivity and regioselectivity predictions, and catalyst design. The process of finding transition states in practice is, however, a laborious multistep operation that requires significant user involvement. Here, we report a highly automated workflow designed to locate transition states for a given elementary reaction with minimal setup overhead. The only essential inputs required from the user are the structures of the separated reactants and products. The seamless workflow combining computational technologies from the fields of cheminformatics, molecular mechanics, and quantum chemistry automatically finds the most probable correspondence between the atoms in the reactants and the products, generates a transition state guess, launches a transition state search through a combined approach involving the relaxing string method and the quadratic synchronous transit, and finally validates the transition state via the analysis of the reactive chemical bonds and imaginary vibrational frequencies as well as by the intrinsic reaction coordinate method. Our approach does not target any specific reaction type, nor does it depend on training data; instead, it is meant to be of general applicability for a wide variety of reaction types. The workflow is highly flexible, permitting modifications such as a choice of accuracy, level of theory, basis set, or solvation treatment. Successfully located transition states can be used for setting up transition state guesses in related reactions, saving computational time and increasing the probability of success. The utility and performance of the method are demonstrated in applications to transition state searches in reactions typical for organic chemistry, medicinal chemistry, and homogeneous catalysis research. In particular, applications of our code to Michael additions, hydrogen abstractions, Diels-Alder cycloadditions, carbene insertions, and an enzyme reaction model involving a molybdenum complex are shown and discussed.

Free Energy Perturbation Calculations of the Thermodynamics of Protein Side-Chain Mutations.

Protein side-chain mutation is fundamental both to natural evolutionary processes and to the engineering of protein therapeutics, which constitute an increasing fraction of important medications. Molecular simulation enables the prediction of the effects of mutation on properties such as binding affinity, secondary and tertiary structure, conformational dynamics, and thermal stability. A number of widely differing approaches have been applied to these predictions, including sequence-based algorithms, knowledge-based potential functions, and all-atom molecular mechanics calculations. Free energy perturbation theory, employing all-atom and explicit-solvent molecular dynamics simulations, is a rigorous physics-based approach for calculating thermodynamic effects of, for example, protein side-chain mutations. Over the past several years, we have initiated an investigation of the ability of our most recent free energy perturbation methodology to model the thermodynamics of protein mutation for two specific problems: protein-protein binding affinities and protein thermal stability. We highlight recent advances in the field and outline current and future challenges.

Acylguanidine Beta Secretase 1 Inhibitors: A Combined Experimental and Free Energy Perturbation Study.

A series of acylguanidine beta secretase 1 (BACE1) inhibitors with modified scaffold and P3 pocket substituent was synthesized and studied with free energy perturbation (FEP) calculations. The resulting molecules showed potencies in enzymatic BACE1 inhibition assays up to 1 nM. The correlation between the predicted activity from the FEP calculations and the experimental activity was good for the P3 pocket substituents. The average mean unsigned error (MUE) between prediction and experiment was 0.68 ± 0.17 kcal/mol for the default 5 ns lambda window simulation time improving to 0.35 ± 0.13 kcal/mol for 40 ns. FEP calculations for the P2' pocket substituents on the same acylguanidine scaffold also showed good agreement with experiment and the results remained stable with repeated simulations and increased simulation time. It proved more difficult to use FEP calculations to study the scaffold modification from increasing 5 to 6 and 7 membered-rings. Although prediction and experiment were in agreement for short 2 ns simulations, as the simulation time increased the results diverged. This was improved by the use of a newly developed "Core Hopping FEP+" approach, which also showed improved stability in repeat calculations. The origins of these differences along with the value of repeat and longer simulation times are discussed. This work provides a further example of the use of FEP as a computational tool for molecular design.

Predicting the Effect of Amino Acid Single-Point Mutations on Protein Stability-Large-Scale Validation of MD-Based Relative Free Energy Calculations.

The stability of folded proteins is critical to their biological function and for the efficacy of protein therapeutics. Predicting the energetic effects of protein mutations can improve our fundamental understanding of structural biology, the molecular basis of diseases, and possible routes to addressing those diseases with biological drugs. Identifying the effect of single amino acid point mutations on the thermodynamic equilibrium between the folded and unfolded states of a protein can pinpoint residues of critical importance that should be avoided in the process of improving other properties (affinity, solubility, viscosity, etc.) and suggest changes at other positions for increasing stability in protein engineering. Multiple computational tools have been developed for in silico predictions of protein stability in recent years, ranging from sequence-based empirical approaches to rigorous physics-based free energy methods. In this work, we show that FEP+, which is a free energy perturbation method based on all-atom molecular dynamics simulations, can provide accurate thermal stability predictions for a wide range of biologically relevant systems. Significantly, the FEP+ approach, while originally developed for relative binding free energies of small molecules to proteins and not specifically fitted for protein stability calculations, performs well compared to other methods that were fitted specifically to predict protein stability. Here, we present the broadest validation of a rigorous free energy-based approach applied to protein stability reported to date: 700+ single-point mutations spanning 10 different protein targets. Across the entire data set, we correctly classify the mutations as stabilizing or destabilizing in 84% of the cases, and obtain statistically significant predictions as compared with experiment [average error of ~1.6kcal/mol and coefficient of determination (R2) of 0.40]. This study demonstrates, for the first time in a large-scale validation, that rigorous free energy calculations can be used to predict changes in protein stability from point mutations without parameterization or system-specific customization, although further improvements should be possible with additional sampling and a better representation of the unfolded state of the protein. Here, we describe the FEP+ method as applied to protein stability calculations, summarize the large-scale retrospective validation results, and discuss limitations of the method, along with future directions for further improvements.

A Comprehensive Computational Study of the Interaction between Human Serum Albumin and Fullerenes.

Human serum albumin (HSA) is the most abundant blood plasma protein, which transports fatty acids, hormones, and drugs. We consider nanoparticle-HSA interactions by investigating the binding of HSA with three fullerene analogs. Long MD simulations, quantum mechanical (fragment molecular orbital, energy decomposition analysis, atoms-in-molecules), and free energy methods elucidated the binding mechanism in these complexes. Such a systematic study is valuable due to the lack of comprehensive theoretical approaches to date. The main elements of the mechanism include the following: binding to IIA site results in allosteric modulation of the IIIA and heme binding sites with an increase in α-helical structure of IIIA. Fullerenes displayed high binding affinities for HSA; therefore, HSA can be used as a fullerene carrier, facilitating any toxic function the fullerene may exert. Complex formation is driven by hydrogen bonding, van der Waals, nonpolar, charge transfer, and dispersion energy contributions. Proper functionalization of C60 has enhanced its binding to HSA by more than an order of magnitude. This feature may be important for biological applications (e.g., photodynamic therapy of cancer). Satisfactory agreement with relevant experimental and theoretical data has been obtained.

Accurate Binding Free Energy Predictions in Fragment Optimization.

Predicting protein-ligand binding free energies is a central aim of computational structure-based drug design (SBDD)--improved accuracy in binding free energy predictions could significantly reduce costs and accelerate project timelines in lead discovery and optimization. The recent development and validation of advanced free energy calculation methods represents a major step toward this goal. Accurately predicting the relative binding free energy changes of modifications to ligands is especially valuable in the field of fragment-based drug design, since fragment screens tend to deliver initial hits of low binding affinity that require multiple rounds of synthesis to gain the requisite potency for a project. In this study, we show that a free energy perturbation protocol, FEP+, which was previously validated on drug-like lead compounds, is suitable for the calculation of relative binding strengths of fragment-sized compounds as well. We study several pharmaceutically relevant targets with a total of more than 90 fragments and find that the FEP+ methodology, which uses explicit solvent molecular dynamics and physics-based scoring with no parameters adjusted, can accurately predict relative fragment binding affinities. The calculations afford R(2)-values on average greater than 0.5 compared to experimental data and RMS errors of ca. 1.1 kcal/mol overall, demonstrating significant improvements over the docking and MM-GBSA methods tested in this work and indicating that FEP+ has the requisite predictive power to impact fragment-based affinity optimization projects.

Stereoselective synthesis and molecular modeling of chiral cyclopentanes.

The reaction of 3-methyseleno-2-methylselenomethyl-propene with benzyl 2,3-anhydro-4-O-triflyl-ß-L-ribopyranoside provides a major convenient enantiomeric product of 1-methylene-(benzyl3,4-dideoxy-α-D-arabinopyranoso)-[3,4-c]-cyclopentane, with benzyl-2,3-anhydro-4-deoxy-4-C-(2-methyl- propen-3-yl)-α-D-lyxopyranoside as a minor product. While the reaction of 3-methyseleno-2-[methylselenomethyl]-propene with benzyl 2,3-anhydro-4-O-triflyl-α-D-ribopyranoside produces a good yield of benzyl-2,3-anhydro-4-deoxy-4-C-(2-methylpropen-3-yl)-α-D-lyxo-pyranoside. Molecular modeling and molecular dynamics simulations indicate that the intermediate in the reaction of the ß-L sugar frequently occupies an optimal conformation that leads to the formation of cyclopentane, while the intermediate in the reaction of the α-D sugar has a very small probability. The results point to the dominant role of the ß-L sugar intermediate in controlling the cyclopentane formation.

Pyrido[4,3-e][1,2,4]triazolo[4,3-a]pyrazines as Selective, Brain Penetrant Phosphodiesterase 2 (PDE2) Inhibitors.

A novel series of pyrido[4,3-e][1,2,4]triazolo[4,3-a]pyrazines is reported as potent PDE2/PDE10 inhibitors with drug-like properties. Selectivity for PDE2 was obtained by introducing a linear, lipophilic moiety on the meta-position of the phenyl ring pending from the triazole. The SAR and protein flexibility were explored with free energy perturbation calculations. Rat pharmacokinetic data and in vivo receptor occupancy data are given for two representative compounds 6 and 12.

Accurate and reliable prediction of relative ligand binding potency in prospective drug discovery by way of a modern free-energy calculation protocol and force field.

Designing tight-binding ligands is a primary objective of small-molecule drug discovery. Over the past few decades, free-energy calculations have benefited from improved force fields and sampling algorithms, as well as the advent of low-cost parallel computing. However, it has proven to be challenging to reliably achieve the level of accuracy that would be needed to guide lead optimization (∼5× in binding affinity) for a wide range of ligands and protein targets. Not surprisingly, widespread commercial application of free-energy simulations has been limited due to the lack of large-scale validation coupled with the technical challenges traditionally associated with running these types of calculations. Here, we report an approach that achieves an unprecedented level of accuracy across a broad range of target classes and ligands, with retrospective results encompassing 200 ligands and a wide variety of chemical perturbations, many of which involve significant changes in ligand chemical structures. In addition, we have applied the method in prospective drug discovery projects and found a significant improvement in the quality of the compounds synthesized that have been predicted to be potent. Compounds predicted to be potent by this approach have a substantial reduction in false positives relative to compounds synthesized on the basis of other computational or medicinal chemistry approaches. Furthermore, the results are consistent with those obtained from our retrospective studies, demonstrating the robustness and broad range of applicability of this approach, which can be used to drive decisions in lead optimization.

3D hydrophobic moment vectors as a tool to characterize the surface polarity of amphiphilic peptides.

The interaction of membranes with peptides and proteins is largely determined by their amphiphilic character. Hydrophobic moments of helical segments are commonly derived from their two-dimensional helical wheel projections, and the same is true for ß-sheets. However, to the best of our knowledge, there exists no method to describe structures in three dimensions or molecules with irregular shape. Here, we define the hydrophobic moment of a molecule as a vector in three dimensions by evaluating the surface distribution of all hydrophilic and lipophilic regions over any given shape. The electrostatic potential on the molecular surface is calculated based on the atomic point charges. The resulting hydrophobic moment vector is specific for the instantaneous conformation, and it takes into account all structural characteristics of the molecule, e.g., partial unfolding, bending, and side-chain torsion angles. Extended all-atom molecular dynamics simulations are then used to calculate the equilibrium hydrophobic moments for two antimicrobial peptides, gramicidin S and PGLa, under different conditions. We show that their effective hydrophobic moment vectors reflect the distribution of polar and nonpolar patches on the molecular surface and the calculated electrostatic surface potential. A comparison of simulations in solution and in lipid membranes shows how the peptides undergo internal conformational rearrangement upon binding to the bilayer surface. A good correlation with solid-state NMR data indicates that the hydrophobic moment vector can be used to predict the membrane binding geometry of peptides. This method is available as a web application on http://www.ibg.kit.edu/HM/.

Ligand photo-isomerization triggers conformational changes in iGluR2 ligand binding domain.

Neurological glutamate receptors bind a variety of artificial ligands, both agonistic and antagonistic, in addition to glutamate. Studying their small molecule binding properties increases our understanding of the central nervous system and a variety of associated pathologies. The large, oligomeric multidomain membrane protein contains a large and flexible ligand binding domains which undergoes large conformational changes upon binding different ligands. A recent application of glutamate receptors is their activation or inhibition via photo-switchable ligands, making them key systems in the emerging field of optochemical genetics. In this work, we present a theoretical study on the binding mode and complex stability of a novel photo-switchable ligand, ATA-3, which reversibly binds to glutamate receptors ligand binding domains (LBDs). We propose two possible binding modes for this ligand based on flexible ligand docking calculations and show one of them to be analogues to the binding mode of a similar ligand, 2-BnTetAMPA. In long MD simulations, it was observed that transitions between both binding poses involve breaking and reforming the T686-E402 protein hydrogen bond. Simulating the ligand photo-isomerization process shows that the two possible configurations of the ligand azo-group have markedly different complex stabilities and equilibrium binding modes. A strong but slow protein response is observed after ligand configuration changes. This provides a microscopic foundation for the observed difference in ligand activity upon light-switching.