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 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.

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.

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.

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.

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.

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.

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/.

Ground state spectroscopy of hydroxyquinolines: evidence for the formation of protonated species in water-rich dioxane-water mixtures.

We have recently used 6-, 7-, and 8-hydroxyquinolines (HQs) as fluorescent probes to study the binding mechanism in one of the drug binding sites of human serum albumin. In the present work we study the absorption spectra of the HQ molecules in neat and binary mixtures of dioxane and water in order to identify the different tautomeric species in the ground state. This study should help in identifying the environment in nanocavities of macromolecules when HQs are used as local reporters. The enol form is shown to be the only tautomer for the three HQs in dioxane and water, with the exception of 7HQ in which both the enol and the zwitterion tautomers exist in equilibrium in water. The results are confirmed by the density functional theory (DFT) calculations using the B3LYP method with a 6-311++G(2d,p) basis set. In water-rich dioxane mixtures, all HQs are protonated. The results were confirmed by comparing the absorption spectra in binary solvents with those in acidic and basic aqueous solutions, and by DFT calculations of the Franck-Condon S1 ← S0 transitions. The number of water molecules solvating the polar sites in each HQ molecule was estimated from the spectral change in the binary solvent mixtures, and structures were calculated by DFT. Mapping the water density around the polar sites in each HQ using molecular dynamics (MD) simulations shows well-defined hydrogen bonds around the N-heteroatom in each HQ molecule. Water density is only well-defined around the hydroxyl group in 8HQ. The MD simulations indicate free rotation of the OH group in 6HQ and 7HQ, and the stability of the cis-isomer in 8HQ. The results point to the unique spectral signatures of 7HQ in water which make this molecule a potential probe to detect the presence of water in nanocavities of macromolecules, and to the ability of the three HQs to detect acidic media in binding sites.

Controlling biological activity with light: diarylethene-containing cyclic peptidomimetics.

Photobiological processes in nature are usually triggered by nonpeptidic chromophores or by modified side chains. A system is presented in which the polypeptide backbone itself can be conformationally switched by light. An amino acid analogue was designed and synthesized based on a reversibly photoisomerizable diarylethene scaffold. This analogue was incorporated into the cyclic backbone of the antimicrobial peptide gramicidin S at several sites. The biological activity of the resulting peptidomimetics could then be effectively controlled by ultraviolet/visible light within strictly defined spatial and temporal limits.

Recent progress in biological charge transfer: theory and simulation.

In this contribution, we discuss three recent developments in atomistic biological charge transfer theory. First, in the context of Marcus' classical theory of charge transfer, key quantities of the theory such as driving forces and reorganization enthalpies are now accessible by thermodynamic integration schemes within standard molecular dynamics simulations at high accuracy. Second, direct simulations of charge transfer enable the computation of fast charge transfer reaction rates without having to resort to Marcus' theory. Finally, exploring the electronic structure beyond that of hitherto presumed centers of localization helps to identify new stepping stones of charge transfer reactions. This article is part of a Special Issue entitled: 17th European Bioenergetics Conference (EBEC 2012).

Rotamer decomposition and protein dynamics: efficiently analyzing dihedral populations from molecular dynamics.

Molecular mechanics methods have matured into powerful methods to understand the dynamics and flexibility of macromolecules and especially proteins. As multinanosecond to microsecond length molecular dynamics (MD) simulations become commonplace, advanced analysis tools are required to generate scientifically useful information from large amounts of data. Some of the key degrees of freedom to understand protein flexibility and dynamics are the amino acid residue side chain dihedral angles. In this work, we present an easily automated way to summarize and understand the relevant dihedral populations. A tremendous reduction in complexity is achieved by describing dihedral timeseries in terms of histograms decomposed into Gaussians. Using the familiar and widely studied protein lysozyme, it is demonstrated that our approach captures essential properties of protein structure and dynamics. A simple classification scheme is proposed that indicates the rotational state population for each dihedral angle of interest and allows a decision if a given side chain or peptide backbone fragment remains rigid during the course of an MD simulation, adopts a converged distribution between conformational substates or has not reached convergence yet.

A contribution to the drug resistance mechanism of darunavir, amprenavir, indinavir, and saquinavir complexes with HIV-1 protease due to flap mutation I50V: a systematic MM-PBSA and thermodynamic integration study.

The emergence of HIV-1 drug-resistant mutations is the major problem against AIDS treatment. We employed molecular dynamics (MD) calculations and free energy (MM-PBSA and thermodynamic integration) analyses on wild-type (WT) and mutated HIV-1 protease (HIV-1 PR) complexes with darunavir, amprenavir, indinavir, and saquinavir to clarify the mechanism of resistance due to the I50V flap mutation. Conformational analysis showed that the protease flaps are increasingly flexible in the I50V complexes. In the WT, stabilization of the HIV-1 PR structure is achieved via interflap and water-mediated hydrogen-bonding interactions between the flaps. Furthermore, hydrogen bonds between drugs and binding cavity residues (Asp29/29'/30/30') are crucial for effective inhibition. All these interactions were significantly diminished (or absent) in the mutated forms, thus denoting their importance toward binding. Thermodynamic integration calculations reproduced the experimental data to within ≈1 kcal mol⁻¹ and showed that the I50V mutation results in weaker binding free energies for all analyzed complexes with respect to the WT. It was observed that the loss in binding energy upon mutation was mostly enthalpically driven in all complexes, with the greatest effect coming from the reduction of van der Waals interactions. Our results motivated us to test two novel compounds that have been synthesized to maximize interactions with HIV-1 PR. MM-PBSA and TI calculations showed that compound 3c (Ghosh et al. Bioorg. Med. Chem. Lett. 2012, 22, 2308) is a promising protease inhibitor, which presents very effective binding to the WT PR (ΔG(MM-PBSA) = -17.2 kcal mol⁻¹, ΔG(exp) = -16.1 kcal mol⁻¹). Upon I50V mutation, the complex binding free energy was weakened by a ΔΔG(TI) of 1.8 kcal mol⁻¹, comparable to the marketed inhibitors. This predicts that I50V may confer low resistance to 3c. This computational comparative study contributes toward elucidation of the I50V drug-resistance mechanism in HIV-1 PR.

Peptide-lipid interactions of the stress-response peptide TisB that induces bacterial persistence.

The bacterial stress-response peptide TisB in Escherichia coli has been suggested to dissipate the transmembrane potential, such that the depletion of ATP levels induces the formation of dormant persister cells which can eventually form biofilms. We studied the structure and membrane interactions of TisB to find out whether it forms pores or other proton-selective channels. Circular dichroism revealed an amphiphilic α-helical structure when reconstituted in lipid vesicles, and oriented circular dichroism showed that the helix assumes a transmembrane alignment. The addition of TisB to dye-loaded vesicles caused leakage only at very high peptide concentration, notably with a Hill coefficient of 2, which suggests that dimers must be involved. Coarse-grained molecular dynamics simulations showed that membrane binding of monomeric TisB is rapid and spontaneous, and transmembrane insertion is energetically feasible. When TisB oligomers are assembled as transmembrane pores, these channels collapse during the simulations, but transmembrane dimers are found to be stable. Given the pattern of charges on the amphiphilic TisB helix, we postulate that antiparallel dimers could be assembled via a ladder of salt bridges. This electrostatic charge-zipper could enable protons to pass along a wire of trapped water molecules across the hydrophobic membrane.

Free energy calculations on the binding of novel thiolactomycin derivatives to E. coli fatty acid synthase I.

Finding novel antibiotics to combat the rise of drug resistance in harmful bacteria is of enormous importance for human health. Computational drug design can be employed to aid synthetic chemists in the search for new potent inhibitors. In recent years, molecular dynamics based free energy calculations have emerged as a useful tool to accurately calculate receptor binding affinities of novel or modified ligands. While being significantly more demanding in computational resources than simpler docking algorithms, they can be employed to obtain reliable estimates of the effect individual functional groups have on protein-ligand complex binding constants. Beta-ketoacyl [acyl carrier protein] synthase I, KAS I, facilitates a critical chain elongation step in the fatty acid synthesis pathway. Since the bacterial type II lipid synthesis system is fundamentally different from the mammalian type I multi-enzyme complex, this enzyme represents a promising target for the design of specific antibiotics. In this work, we study the binding of several recently synthesized derivatives of the natural KAS I inhibitor thiolactomycin in detail based on atomistic modeling. From extensive thermodynamic integration calculations the effect of changing functional groups on the thiolactone scaffold was determined. Four ligand modifications were predicted to show improved binding to the E. coli enzyme, pointing the way towards the design of thiolactomycin derivatives with binding constants in the nanomolar range.

Molecular origin of differences in hole and electron mobility in amorphous Alq3--a multiscale simulation study.

In order to determine the molecular origin of the difference in electron and hole mobilities of amorphous thin films of Alq(3) (meridional Alq(3) (tris(8-hydroxyquinoline) aluminium)) we performed multiscale simulations covering quantum mechanics, molecular mechanics and lattice models. The study includes realistic disordered morphologies, polarized site energies to describe diagonal disorder, quantum chemically calculated transfer integrals for the off-diagonal disorder, inner sphere reorganization energies and an approximative scheme for outer sphere reorganization energies. Intermolecular transfer rates were calculated via Marcus-theory and mobilities were simulated via kinetic Monte Carlo simulations and by a Master Equation approach. The difference in electron and hole mobility originates from the different localization of charge density in the radical anion (more delocalized) compared to the radical cation (more confined). This results in higher diagonal disorder for holes and less favourable overlap properties for the hole transfer integrals leading to an overall higher electron mobility.

Soft-core potentials in thermodynamic integration: comparing one- and two-step transformations.

Molecular dynamics-based free energy calculations allow the determination of a variety of thermodynamic quantities from computer simulations of small molecules. Thermodynamic integration (TI) calculations can suffer from instabilities during the creation or annihilation of particles. This "singularity" problem can be addressed with "soft-core" potential functions which keep pairwise interaction energies finite for all configurations and provide smooth free energy curves. "One-step" transformations, in which electrostatic and van der Waals forces are simultaneously modified, can be simpler and less expensive than "two-step" transformations in which these properties are changed in separate calculations. Here, we study solvation free energies for molecules of different hydrophobicity using both models. We provide recommended values for the two parameters α(LJ) and ß(C) controlling the behavior of the soft-core Lennard-Jones and Coulomb potentials and compare one- and two-step transformations with regard to their suitability for numerical integration. For many types of transformations, the one-step procedure offers a convenient and accurate approach to free energy estimates.

Tautomerism in 7-hydroxyquinoline: a combined experimental and theoretical study in water.

Tautomerism in the ground and excited states of 7-hydroxyquinoline (7HQ) was studied in different solvents using steady-state and lifetime spectroscopic measurements, density functional theory (DFT) calculations, and molecular dynamics (MD) simulations. Equilibrium between the enol and the keto/zwitterion tautomers exists in 7HQ, which is solvent-dependent. Of the solvents used in this study, only in water does the absorbance spectrum of 7HQ show absorption from both the enol and zwitterion tautomers. In addition, in aqueous media, fluorescence is observed from the zwitterion tautomer only, which is attributed to self-quenching of the enol fluorescence by energy transfer to the ground-state zwitterion tautomer and energetically favorable excited-state proton transfer. Solvation of the hydrogen bonding sites of 7HQ was studied in binary mixtures of 1,4-dioxane and water, and three water molecules were estimated to connect the polar sites and induce intermolecular proton transfer. The results are confirmed by DFT calculations showing that three water molecules are the minimum number required to form a stable solvent wire. Mapping the water density around the polar sites using MD simulations shows well-defined hydrogen bonds around the amino and hydroxyl groups of the enol tautomer and slightly less well-defined hydrogen bonds for the zwitterion tautomer. The presence of three-member water wires connecting the polar centers in 7HQ is evident in the MD simulations. The results point to the unique spectral signatures of 7HQ in water, which make this molecule a potential probe to detect the presence of water in nanocavities of macromolecules.

The thermodynamics of charge transfer in DNA photolyase: using thermodynamic integration calculations to analyse the kinetics of electron transfer reactions.

DNA Photolyases are light sensitive oxidoreductases present in many organisms that participate in the repair of photodamaged DNA. They are capable of electron transfer between a bound cofactor and a chain of tryptophan amino acid residues. Due to their unique mechanism and important function, photolyases have been subject to intense study in recent times, with both experimental and computational efforts. In this work, we present a novel application of classical molecular dynamics based free energy calculations, combined with quantum mechanical computations, to biomolecular charge transfer. Our approach allows for the determination of all reaction parameters in Marcus' theory of charge transport. We were able to calculate the free energy profile for the movement of a positive charge along protein sidechains involved in the biomolecule's function as well as charge-transfer rates that are in good agreement with experimental results. Our approach to simulate charge-transfer reactions explicitly includes the influence of protein flexibility and solvent dynamics on charge-transfer energetics. As applied here to a biomolecular system of considerable scientific interest, we believe the method to be easily adaptable to the study of charge-transfer phenomena in biochemistry and other fields.

Direct simulation of electron transfer reactions in DNA radical cations.

The electron transfer properties of DNA radical cations are important in DNA damage and repair processes. Fast long-range charge transfer has been demonstrated experimentally, but the subtle influences that experimental conditions as well as DNA sequences and geometries have on the details of electron transfer parameters are still poorly understood. In this work, we employ an atomistic QM/MM approach, based on a one-electron tight binding Hamiltonian and a classical molecular mechanics forcefield, to conduct nanosecond length MD simulations of electron holes in DNA oligomers. Multiple spontaneous electron transfer events were observed in 100 ns simulations with neighboring adenine or guanine bases. Marcus parameters of charge transfer could be extracted directly from the simulations. The reorganization energy lambda for hopping between neighboring bases was found to be ca. 25 kcal/mol and charge transfer rates of 4.1 x 10(9) s(-1) for AA hopping and 1.3 x 10(9) s(-1) for GG hopping were obtained.