Reweighting from Molecular Mechanics Force Fields to the ANI-2x Neural Network Potential.

To achieve chemical accuracy in free energy calculations, it is necessary to accurately describe the system's potential energy surface and efficiently sample configurations from its Boltzmann distribution. While neural network potentials (NNPs) have shown significantly higher accuracy than classical molecular mechanics (MM) force fields, they have a limited range of applicability and are considerably slower than MM potentials, often by orders of magnitude. To address this challenge, Rufa et al. [Rufa et al. bioRxiv 2020, 10.1101/2020.07.29.227959.] suggested a two-stage approach that uses a fast and established MM alchemical energy protocol, followed by reweighting the results using NNPs, known as endstate correction or indirect free energy calculation. This study systematically investigates the accuracy and robustness of reweighting from an MM reference to a neural network target potential (ANI-2x) for an established data set in vacuum, using single-step free-energy perturbation (FEP) and nonequilibrium (NEQ) switching simulation. We assess the influence of longer switching lengths and the impact of slow degrees of freedom on outliers in the work distribution and compare the results to those of multistate equilibrium free energy simulations. Our results demonstrate that free energy calculations between NNPs and MM potentials should be preferably performed using NEQ switching simulations to obtain accurate free energy estimates. NEQ switching simulations between the MM potentials and NNPs are efficient, robust, and trivial to implement.

Calculations of Absolute Solvation Free Energies with Transformato─Application to the FreeSolv Database Using the CGenFF Force Field.

We recently introduced transformato, an open-source Python package for the automated setup of large-scale calculations of relative solvation and binding free energy differences. Here, we extend the capabilities of transformato to the calculation of absolute solvation free energy differences. After careful validation against the literature results and reference calculations with the PERT module of CHARMM, we used transformato to compute absolute solvation free energies for most molecules in the FreeSolv database (621 out of 642). The force field parameters were obtained with the program cgenff (v2.5.1), which derives missing parameters from the CHARMM general force field (CGenFF v4.6). A long-range correction for the Lennard-Jones interactions was added to all computed solvation free energies. The mean absolute error compared to the experimental data is 1.12 kcal/mol. Our results allow a detailed comparison between the AMBER and CHARMM general force fields and provide a more in-depth understanding of the capabilities and limitations of the CGenFF small molecule parameters.

Protex-A Python utility for proton exchange in molecular dynamics simulations.

Protex is an open-source program that enables proton exchanges of solvent molecules during molecular dynamics simulations. While conventional molecular dynamics simulations do not allow for bond breaking or formation, protex offers an easy-to-use interface to augment these simulations and define multiple proton sites for (de-)protonation using a single topology approach with two different λ-states. Protex was successfully applied to a protic ionic liquid system, where each molecule is prone to (de-)protonation. Transport properties were calculated and compared to experimental values and simulations without proton exchange.

POLθ processes ssDNA gaps and promotes replication fork progression in BRCA1-deficient cells.

Polymerase theta (POLθ) is an error-prone DNA polymerase whose loss is synthetically lethal in cancer cells bearing breast cancer susceptibility proteins 1 and 2 (BRCA1/2) mutations. To investigate the basis of this genetic interaction, we utilized a small-molecule inhibitor targeting the POLθ polymerase domain. We found that POLθ processes single-stranded DNA (ssDNA) gaps that emerge in the absence of BRCA1, thus promoting unperturbed replication fork progression and survival of BRCA1 mutant cells. A genome-scale CRISPR-Cas9 knockout screen uncovered suppressors of the functional interaction between POLθ and BRCA1, including NBN, a component of the MRN complex, and cell-cycle regulators such as CDK6. While the MRN complex nucleolytically processes ssDNA gaps, CDK6 promotes cell-cycle progression, thereby exacerbating replication stress, a feature of BRCA1-deficient cells that lack POLθ activity. Thus, ssDNA gap formation, modulated by cell-cycle regulators and MRN complex activity, underlies the synthetic lethality between POLθ and BRCA1, an important insight for clinical trials with POLθ inhibitors.

Relative binding free energy calculations with transformato: A molecular dynamics engine-independent tool.

We present the software package transformato for the setup of large-scale relative binding free energy calculations. Transformato is written in Python as an open source project (https://github.com/wiederm/transformato); in contrast to comparable tools, it is not closely tied to a particular molecular dynamics engine to carry out the underlying simulations. Instead of alchemically transforming a ligand L 1 directly into another L 2, the two ligands are mutated to a common core. Thus, while dummy atoms are required at intermediate states, in particular at the common core state, none are present at the physical endstates. To validate the method, we calculated 76 relative binding free energy differences Δ Δ G L 1 → L 2 b i n d for five protein-ligand systems. The overall root mean squared error to experimental binding free energies is 1.17 kcal/mol with a Pearson correlation coefficient of 0.73. For selected cases, we checked that the relative binding free energy differences between pairs of ligands do not depend on the choice of the intermediate common core structure. Additionally, we report results with and without hydrogen mass reweighting. The code currently supports OpenMM, CHARMM, and CHARMM/OpenMM directly. Since the program logic to choose and construct alchemical transformation paths is separated from the generation of input and topology/parameter files, extending transformato to support additional molecular dynamics engines is straightforward.

Improving Small Molecule pK a Prediction Using Transfer Learning With Graph Neural Networks.

Enumerating protonation states and calculating microstate pK a values of small molecules is an important yet challenging task for lead optimization and molecular modeling. Commercial and non-commercial solutions have notable limitations such as restrictive and expensive licenses, high CPU/GPU hour requirements, or the need for expert knowledge to set up and use. We present a graph neural network model that is trained on 714,906 calculated microstate pK a predictions from molecules obtained from the ChEMBL database. The model is fine-tuned on a set of 5,994 experimental pK a values significantly improving its performance on two challenging test sets. Combining the graph neural network model with Dimorphite-DL, an open-source program for enumerating ionization states, we have developed the open-source Python package pkasolver, which is able to generate and enumerate protonation states and calculate pK a values with high accuracy.

Alchemical free energy simulations without speed limits. A generic framework to calculate free energy differences independent of the underlying molecular dynamics program.

We describe the theory of the so-called common-core/serial-atom-insertion (CC/SAI) approach to compute alchemical free energy differences and its practical implementation in a Python package called Transformato. CC/SAI is not tied to a specific biomolecular simulation program and does not rely on special purpose code for alchemical transformations. To calculate the alchemical free energy difference between several small molecules, the physical end-states are mutated into a suitable common core. Since this only requires turning off interactions, the setup of intermediate states is straightforward to automate. Transformato currently supports CHARMM and OpenMM as back ends to carry out the necessary molecular dynamics simulations, as well as post-processing calculations. We validate the method by computing a series of relative solvation free energy differences.

Synthesis and dopamine receptor binding of dihydrexidine and SKF 38393 catecholamine-based analogues.

Dopamine is an important neurotransmitter that regulates numerous essential functions, including cognition and voluntary movement. As such, it serves as an important scaffold for synthesis of novel analogues as part of drug development effort to obtain drugs for treatment of neurodegenerative diseases, such as Parkinson's disease. To that end, similarity search of the ZINC database based on two known dopamine-1 receptor (D1R) agonists, dihydrexidine (DHX) and SKF 38393, respectively, was used to predict novel chemical entities with potential binding to D1R. Three compounds that showed the highest similarity index were selected for synthesis and bioactivity profiling. All main synthesis products as well as the isolated intermediates, were properly characterized. The physico-chemical analyses were performed using HRESIMS, GC/MS, LC/MS with UV-Vis detection, and FTIR, 1H NMR and 13C NMR spectroscopy. Binding to D1 and D2 receptors and inhibition of dopamine reuptake via dopamine transporter were measured for the synthesized analogues of DHX and SKF 38393.

Improved Lipophilicity and Aqueous Solubility Prediction with Composite Graph Neural Networks.

The accurate prediction of molecular properties, such as lipophilicity and aqueous solubility, are of great importance and pose challenges in several stages of the drug discovery pipeline. Machine learning methods, such as graph-based neural networks (GNNs), have shown exceptionally good performance in predicting these properties. In this work, we introduce a novel GNN architecture, called directed edge graph isomorphism network (D-GIN). It is composed of two distinct sub-architectures (D-MPNN, GIN) and achieves an improvement in accuracy over its sub-architectures employing various learning, and featurization strategies. We argue that combining models with different key aspects help make graph neural networks deeper and simultaneously increase their predictive power. Furthermore, we address current limitations in assessment of deep-learning models, namely, comparison of single training run performance metrics, and offer a more robust solution.

Fitting quantum machine learning potentials to experimental free energy data: predicting tautomer ratios in solution.

The computation of tautomer ratios of druglike molecules is enormously important in computer-aided drug discovery, as over a quarter of all approved drugs can populate multiple tautomeric species in solution. Unfortunately, accurate calculations of aqueous tautomer ratios-the degree to which these species must be penalized in order to correctly account for tautomers in modeling binding for computer-aided drug discovery-is surprisingly difficult. While quantum chemical approaches to computing aqueous tautomer ratios using continuum solvent models and rigid-rotor harmonic-oscillator thermochemistry are currently state of the art, these methods are still surprisingly inaccurate despite their enormous computational expense. Here, we show that a major source of this inaccuracy lies in the breakdown of the standard approach to accounting for quantum chemical thermochemistry using rigid rotor harmonic oscillator (RRHO) approximations, which are frustrated by the complex conformational landscape introduced by the migration of double bonds, creation of stereocenters, and introduction of multiple conformations separated by low energetic barriers induced by migration of a single proton. Using quantum machine learning (QML) methods that allow us to compute potential energies with quantum chemical accuracy at a fraction of the cost, we show how rigorous relative alchemical free energy calculations can be used to compute tautomer ratios in vacuum free from the limitations introduced by RRHO approximations. Furthermore, since the parameters of QML methods are tunable, we show how we can train these models to correct limitations in the underlying learned quantum chemical potential energy surface using free energies, enabling these methods to learn to generalize tautomer free energies across a broader range of predictions.

Dummy Atoms in Alchemical Free Energy Calculations.

In calculations of relative free energy differences, the number of atoms of the initial and final states is rarely the same. This necessitates the introduction of dummy atoms. These placeholders interact with the physical system only by bonded energy terms. We investigate the conditions necessary so that the presence of dummy atoms does not influence the result of a relative free energy calculation. On the one hand, one has to ensure that dummy atoms only give a multiplicative contribution to the partition function so that their contribution cancels from double-free energy differences. On the other hand, the bonded terms used to attach a dummy atom (or group of dummy atoms) to the physical system have to maintain it in a well-defined position and orientation relative to the physical system. A detailed theoretical analysis of both aspects is provided, illustrated by 24 calculations of relative solvation free energy differences, for which all four legs of the underlying thermodynamic cycle were computed. Cycle closure (or lack thereof) was used as a sensitive indicator to probing the effects of dummy atom treatment on the resulting free energy differences. We find that a naive (but often practiced) treatment of dummy atoms results in errors of up to kBT when calculating the relative solvation free energy difference between two small solutes, such as methane and ammonia. While our analysis focuses on the so-called single topology approach to set up alchemical transformations, similar considerations apply to dual topology, at least many widely used variants thereof.

Structure-Activity Relationships of Novel Thiazole-Based Modafinil Analogues Acting at Monoamine Transporters.

Atypical dopamine reuptake inhibitors, such as modafinil, are used for the treatment of sleeping disorders and investigated as potential therapeutics against cocaine addiction and for cognitive enhancement. Our continuous effort to find modafinil analogues with higher inhibitory activity on and selectivity toward the dopamine transporter (DAT) has previously led to the promising thiazole-containing derivatives CE-103, CE-111, CE-123, and CE-125. Here, we describe the synthesis and activity of a series of compounds based on these scaffolds, which resulted in several new selective DAT inhibitors and gave valuable insights into the structure-activity relationships. Introduction of the second chiral center and subsequent chiral separations provided all four stereoisomers, whereby the S-configuration on both generally exerted the highest activity and selectivity on DAT. The representative compound of this series was further characterized by in silico, in vitro, and in vivo studies that have demonstrated both safety and efficacy profile of this compound class.

Computational Identification of Novel Kir6 Channel Inhibitors.

KATP channels consist of four Kir6.x pore-forming subunits and four regulatory sulfonylurea receptor (SUR) subunits. These channels couple the metabolic state of the cell to membrane excitability and play a key role in physiological processes such as insulin secretion in the pancreas, protection of cardiac muscle during ischemia and hypoxic vasodilation of arterial smooth muscle cells. Abnormal channel function resulting from inherited gain or loss-of-function mutations in either the Kir6.x and/or SUR subunits are associated with severe diseases such as neonatal diabetes, congenital hyperinsulinism, or Cantú syndrome (CS). CS is an ultra-rare genetic autosomal dominant disorder, caused by dominant gain-of-function mutations in SUR2A or Kir6.1 subunits. No specific pharmacotherapeutic treatment options are currently available for CS. Kir6 specific inhibitors could be beneficial for the development of novel drug therapies for CS, particular for mutations, which lack high affinity for sulfonylurea inhibitor glibenclamide. By applying a combination of computational methods including atomistic MD simulations, free energy calculations and pharmacophore modeling, we identified several novel Kir6.1 inhibitors, which might be possible candidates for drug repurposing. The in silico predictions were confirmed using inside/out patch-clamp analysis. Importantly, Cantú mutation C166S in Kir6.2 (equivalent to C176S in Kir6.1) and S1020P in SUR2A, retained high affinity toward the novel inhibitors. Summarizing, the inhibitors identified in this study might provide a starting point toward developing novel therapies for Cantú disease.

Discovery of Potent Inhibitors for the Large Neutral Amino Acid Transporter 1 (LAT1) by Structure-Based Methods.

The large neutral amino acid transporter 1 (LAT1) is a promising anticancer target that is required for the cellular uptake of essential amino acids that serve as building blocks for cancer growth and proliferation. Here, we report a structure-based approach to identify chemically diverse and potent inhibitors of LAT1. First, a homology model of LAT1 that is based on the atomic structures of the prokaryotic homologs was constructed. Molecular docking of nitrogen mustards (NMs) with a wide range of affinity allowed for deriving a common binding mode that could explain the structure-activity relationship pattern in NMs. Subsequently, validated binding hypotheses were subjected to molecular dynamics simulation, which allowed for extracting a set of dynamic pharmacophores. Finally, a library of ~1.1 million molecules was virtually screened against these pharmacophores, followed by docking. Biological testing of the 30 top-ranked hits revealed 13 actives, with the best compound showing an IC50 value in the sub-µM range.

The Use of Dynamic Pharmacophore in Computer-Aided Hit Discovery: A Case Study.

In recent years pharmacophore modeling has become increasingly popular due to the development of software solutions and improvement in algorithms that allowed researchers to focus on interactions between protein and ligands instead of technical details of the software. At the same time, progress in computer hardware made molecular dynamics (MD) simulations on regular PC hardware possible. MD simulations are usually used, within the virtual screening process, to take into account the flexibility of the target and studying it in more realistic way. In order to do so, it is customary to use simulations before the virtual screening process and then use them for collecting some specific conformation of the target used. Furthermore, some researchers have demonstrated that the use of multiple crystal structures of the same protein can be valuable to better explore the role of the ligand within the binding pocket and then evaluate the most important interactions that are created during the host-guest recognition process. Findings derived from the MD analysis, especially focused on interactions, can be in fact exploited as features for pharmacophore generation or constraints to be used in the molecular docking as integrated steps of the whole virtual screening process. In this chapter, we will present the recent advances in the field pharmacophore modeling combined with the use of MD, a field well explored by our research group in the last 2 years.

SAR-Guided Scoring Function and Mutational Validation Reveal the Binding Mode of CGS-8216 at the α1+/γ2- Benzodiazepine Site.

The structural resolution of a bound ligand-receptor complex is a key asset to efficiently drive lead optimization in drug design. However, structural resolution of many drug targets still remains a challenging endeavor. In the absence of structural knowledge, scientists resort to structure-activity relationships (SARs) to promote compound development. In this study, we incorporated ligand-based knowledge to formulate a docking scoring function that evaluates binding poses for their agreement with a known SAR. We showcased this protocol by identifying the binding mode of the pyrazoloquinolinone (PQ) CGS-8216 at the benzodiazepine binding site of the GABAA receptor. Further evaluation of the final pose by molecular dynamics and free energy simulations revealed a close proximity between the pendent phenyl ring of the PQ and γ2D56, congruent with the low potency of carboxyphenyl analogues. Ultimately, we introduced the γ2D56A mutation and in fact observed a 10-fold potency increase in the carboxyphenyl analogue, providing experimental evidence in favor of our binding hypothesis.

A daily single dose of a novel modafinil analogue CE-123 improves memory acquisition and memory retrieval.

Dopamine reuptake inhibitors have been shown to improve cognitive parameters in various tasks and animal models. We recently reported a series of modafinil analogues, of which the most promising, 5-((benzhydrylsulfinyl)methyl) thiazole (CE-123), was selected for further development. The present study aims to characterize pharmacological properties of CE-123 and to investigate the potential to enhance memory performance in a rat model. In vitro transporter assays were performed in cells expressing human transporters. CE-123 blocked uptake of [3H] dopamine (IC50 = 4.606 µM) while effects on serotonin (SERT) and the norepinephrine transporter (NET) were negligible. Blood-brain barrier and pharmacokinetic studies showed that the compound reached the brain and lower elimination than R-modafinil. The Pro-cognitive effect was evaluated in a spatial hole-board task in male Sprague-Dawley rats and CE-123 enhances memory acquisition and memory retrieval, represented by significantly increased reference memory indices and shortened latency. Since DAT blockers can be considered as indirect dopamine receptor agonists, western blotting was used to quantify protein levels of dopamine receptors D1R, D2R and D5R and DAT in the synaptosomal fraction of hippocampal subregions CA1, CA3 and dentate gyrus (DG). CE-123 administration in rats increased total DAT levels and D1R protein levels were significantly increased in CA1 and CA3 in treated/trained groups. The increase of D5R was observed in DG only. Dopamine receptors, particularly D1R, seem to play a role in mediating CE-123-induced memory enhancement. Dopamine reuptake inhibition by CE-123 may represent a novel and improved stimulant therapeutic for impairments of cognitive functions.

Heterocyclic Analogues of Modafinil as Novel, Atypical Dopamine Transporter Inhibitors.

Modafinil is a wake promoting compound with high potential for cognitive enhancement. It is targeting the dopamine transporter (DAT) with moderate selectivity, thereby leading to reuptake inhibition and increased dopamine levels in the synaptic cleft. A series of modafinil analogues have been reported so far, but more target-specific analogues remain to be discovered. It was the aim of this study to synthesize and characterize such analogues and, indeed, a series of compounds were showing higher activities on the DAT and a higher selectivity toward DAT versus serotonin and norepinephrine transporters than modafinil. This was achieved by substituting the amide moiety by five- and six-membered aromatic heterocycles. In vitro studies indicated binding to the cocaine pocket on DAT, although molecular dynamics revealed binding different from that of cocaine. Moreover, no release of dopamine was observed, ruling out amphetamine-like effects. The absence of neurotoxicity of a representative analogue may encourage further preclinical studies of the above-mentioned compounds.

A Molecular Dynamics-Shared Pharmacophore Approach to Boost Early-Enrichment Virtual Screening: A Case Study on Peroxisome Proliferator-Activated Receptor†α.

Molecular dynamics (MD) simulations can be used, prior to virtual screening, to add flexibility to proteins and study them in a dynamic way. Furthermore, the use of multiple crystal structures of the same protein containing different co-crystallized ligands can help elucidate the role of the ligand on a protein's active conformation, and then explore the most common interactions between small molecules and the receptor. In this work, we evaluated the contribution of the combined use of MD on crystal structures containing the same protein but different ligands to examine the crucial ligand-protein interactions within the complexes. The study was carried out on peroxisome proliferator-activated receptor α (PPARα). Findings derived from the dynamic analysis of interactions were then used as features for pharmacophore generation and constraints for generating the docking grid for use in virtual screening. We found that information derived from short multiple MD simulations using different molecules within the binding pocket of the target can improve the early enrichment of active ligands in the virtual screening process for this receptor. In the end we adopted a consensus scoring based on docking score and pharmacophore alignment to rank our dataset. Our results showed an improvement in virtual screening performance in early recognition when screening was performed with the Molecular dYnamics SHAred PharmacophorE (MYSHAPE) approach.

Common Hits Approach: Combining Pharmacophore Modeling and Molecular Dynamics Simulations.

We present a new approach that incorporates flexibility based on extensive MD simulations of protein-ligand complexes into structure-based pharmacophore modeling and virtual screening. The approach uses the multiple coordinate sets saved during the MD simulations and generates for each frame a pharmacophore model. Pharmacophore models with the same pharmacophore features are pooled. In this way the high number of pharmacophore models that results from the MD simulation is reduced to only a few hundred representative pharmacophore models. Virtual screening runs are performed with every representative pharmacophore model; the screening results are combined and rescored to generate a single hit-list. The score for a particular molecule is calculated based on the number of representative pharmacophore models which classified it as active. Hence, the method is called common hits approach (CHA). The steps between the MD simulation and the final hit-list are performed automatically and without user interaction. We test the performance of CHA for virtual screening using screening databases with active and inactive compounds for 40 protein-ligand systems. The results of the CHA are compared to the (i) median screening performance of all representative pharmacophore models of protein-ligand systems, as well as to the virtual screening performance of (ii) a random classifier, (iii) the pharmacophore model derived from the experimental structure in the PDB, and (iv) the representative pharmacophore model appearing most frequently during the MD simulation. For the 34 (out of 40) protein-ligand complexes, for which at least one of the approaches was able to perform better than a random classifier, the highest enrichment was achieved using CHA in 68% of the cases, compared to 12% for the PDB pharmacophore model and 20% for the representative pharmacophore model appearing most frequently. The availabilithy of diverse sets of different pharmacophore models is utilized to analyze some additional questions of interest in 3D pharmacophore-based virtual screening.